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SUMMARY 


This  report  summarizes  research  done  in  the  last  two  and  one-half  years  as  an  employee  of 
the  Operational  Technologies  Corporation.  The  original  assignment  was  with  the  Radiation 
Analysis  Branch  of  the  Directed  Energy  Division  (AL/OEDA)  but,  due  to  a  number  of 
reorganizations,  was  finished  with  the  Biomathematics  Branch  of  the  Mathematical  Products 
Division  (AL/OESB). 

This  effort  addressed  the  direct  problem  in  electromagnetic  wave  propagation  and  scat¬ 
tering  in  dispersive  dielecfrics  objects.  Particularly,  the  research  concentrated  on  the  devel¬ 
opment  of  numerical  approaches  that  can  solve  the  problem  in  a  timely  fashion  as  accurately 
as  possible  given  a  reasonable  amount  of  computational  resources.  In  addition,  analytical 
methods  were  developed  for  the  purpose  of  extracting  features  of  the  problem  which  are 
specicil  to  dispersive  objects. 

First,  four  existing  Finite  Difference  Time  Domain  extensions  for  the  modeling  of  pulse 
propagation  in  Debye  or  Lorentz  dispersive  media  were  analyzed  through  studying  the  stabil¬ 
ity  and  phase  error  properties  of  the  coupled  difference  equations  corresponding  to  Maxwell’s 
equations  and  to  the  equations  for  the  dispersion.  For  good  overall  accuracy  we  showed  that 
all  schemes  should  be  run  at  their  Courant  stability  limit,  and  that  the  timestep  should 
finely  resolve  the  medium  timescales.  Particularly,  for  Debye  schemes  it  should  be  at  least 
At  =  10"^r,  while  for  Lorentz  schemes  it  should  be  At  =  10"V,  where  r  is  a  typical 
medium  relaxation  time.  A  numerical  experiment  with  a  Debye  medium  confirmed  this.  We 
also  determined  that  two  of  the  discretizations  for  Debye  media  are  totally  equivalent.  In  the 
Lorentz  medium  case  we  established  that  the  method  that  uses  the  polarization  differential 
equation  to  model  dispersion  is  stable  for  all  wavenumbers,  and  that  the  method  using  the 
local-in-time  constitutive  relation  is  weakly  unstable  for  modes  with  wavenumber  k  such  that 
fcAx  >  7r/2. 

In  order  to  understand  more  fuUy  the  discretization  requirements  and  the  general  behavior 
of  existing  numerical  methods  for  dispersive  media  we  considered  the  propagation  of  arbitrary 
electromagnetic  pulses  in  anomalously  dispersive  dielectrics  characterized  by  M  relaxation 
processes.  A  partial  differential  equation  for  the  electric  field  in  the  dielectric  was  derived 
and  anzJyzed.  This  single  equation  describes  a  hierarchy  of  M  +  1  wave  types,  each  type 
characterized  by  an  attenuation  coefficient  and  a  wave  speed.  Our  analysis  identified  a  skin- 
depth”  where  the  pulse  response  is  described  by  a  telegrapher’s  equation  with  smoothing 


terms,  travels  with  the  wavefront  speed,  and  decays  exponentially.  Past  this  shallow  depth 
we  showed  that  the  pulse  response  is  described  by  a  weakly  dispersive  advection-diffusion 
equation,  travels  with  the  sub-chajacteristic  advection  speed  equal  to  the  zero-frequency 
phase  velocity  in  the  dielectric,  and  decays  cdgebraically.  The  analysis  was  verified  with  a 
numericaJ  simulation.  The  relevance  of  our  results  to  the  development  of  numericcJ  methods 
for  such  problems  was  discussed. 

Armed  with  the  knowledge  that  finite  difference  schemes,  which  are  more  accurate  in 
space  than  in  time,  are  better  suited  for  modeling  pulse  propagation  in  dispersive  dielectrics 
we  next  determined  some  important  properties  of  such  high-order  schemes  in  comparison  to 
standard  schemes.  For  Finite  Difference  Time  Domain  methods  we  determined  the  spatial 
resolution  of  the  discretized  domain  in  terms  of  the  total  computation  time  and  the  desired 
phase  error.  It  was  shown  that  the  spatial  step  should  vary  as  Ax  ~  g  in  order  to 

maintain  a  prescribed  phase  error  level  throughout  the  computation  time  tcy  where  s  (=2 
or  4)  is  the  spatial  order  of  accuracy  of  the  scheme  and  ^  is  a  geometric  factor.  Significantly, 
we  showed  that  the  thumb  rule  of  using  10-20  points  per  wavelength  to  determine  the  spatial 
cell  size  for  the  standcird  scheme  is  not  optimal.  Our  results  were  verified  by  numericzJ 
simulations  in  two  dimensions  with  the  Yee  scheme  and  the  new  4th-order  accurate  scheme. 

Finally,  we  examined  further  aspects  of  electromagnetic  pulse  propagation  in  anomalously 
dispersive  media,  using  the  Debye  model  as  an  example.  Short-pulse,  long-pulse,  short- 
time,  and  long-time  approximations  and  amplitude  rate  of  decay  estimates  were  derived 
with  asymptotic  methods.  We  cdso  studied  the  following  problem:  Knowing  only  the  peak 
amplitude  cind  power  density  of  an  incident  pulse,  what  can  be  said  about  the  peak  amplitude 
of  the  propagated  pulse?  We  provided  sharp  upper  and  lower  bounds  for  the  propagated 
amplitude  which  may  be  useful  in  controlling  the  electromagnetic  interference  or  the  damage 
produced  in  dispersive  media.  We  explained  a  factor-of-nine  effect  in  the  speed  of  waves 
in  a  Debye  model  for  water,  which  seems  to  have  been  previously  unnoticed,  and  some 
observations  from  experimentad  studies  of  pulse  propagation  in  biological  materials.  FinaJly, 
we  proposed  some  guidelines  for  sample  size  in  Transmission  Time  Domain  Spectroscopy 
studies  of  dielectrics. 


stability  and  Phase  Error  Analysis  of 
FD-TD  in  Dispersive  Dielectrics 


Abstract — Four  FD-TD  extension*  for  the  modeling  of  pulse 
propagation  in  Debye  or  Lorent*  dispersive  media  are  ana¬ 
lysed  through  studying  the  stability  and  phase  error  prop¬ 
erties  of  the  coupled  difference  equations  corresponding  to 
Maxwell’s  equations  and  to  the  equations  for  the  dispersion. 
For  good  overall  accuracy  we  show  tha#  all  schemes  should 
be  run  at  their  Courant  stability  limit,  and  that  the  timestep 
should  finely  resolve  the  medium  timescales.  Pwticularly, 
for  Debye  schemes  It  should  be  at  least  At  =  lo  while 
for  Lorents  schemes  it  should  be  At  =  lo  *t,  where  r  is  a 
typical  medium  relaxation  time.  A  numerical  experiment 
with  a  Debye  medium  confirms  this.  ^Ve  have  determined 
that  two  of  the  discretisations  for  Debye  media  are  totally 
equivalent.  In  the  Lorents  medium  case  we  establish  that 
the  method  that  uses  the  polarisation  differential  equation 
to  model  dispersion  is  stable  for  all  wavenumbers,  and  that 
the  method  using  the  local-in-time  constitutive  relation  is 
weakly  unstable  for  modes  with  wavenumber  k  such  that 
kAx  >  -x/a. 

Keywords —  Dispersive  media,  finite-difference  time-domain 
method,  algorithm  stability  and  phase  error,  RF  dosimetry, 
bioelectromagnetics. 

1.  Introduction 

THE  Computational  Electromagnetics  community  has  de- 
vcloped  a  variety  of  extensions  of  the  popular  FD-TD  nu¬ 
merical  method  [l]-[2]  for  the  modeling  of  pulse  propaga¬ 
tion  in  dispersive  media  with  complex  geometry.  A  repre¬ 
sentative  list  of  these  extensions  is  [3]- [7].  Some  of  these 
extensions  append  to  the  Partial  Differential  Equations 
(PDEs),  that  express  Maxwell’s  equations,  a  set  of  Ordi¬ 
nary  Differential  Equations  (ODEs)  that  either  describe 
the  local-in-time  constitutive  relation  [3]  involving  the  dis¬ 
placement  D  and  the  electric  field  E,  or  the  dynamic  evolu¬ 
tion  of  the  polarization  P  excited  by  the  propagating  elec¬ 
tric  field  [4]-[5].  Other  popular  extensions  [6]-[7]  use  a  con¬ 
volution  representation  of  the  constitutive  relation  which 
is  updated  in  sync  with  the  time  update  of  the  FD-TD  dis¬ 
cretized  PDEs.  In  the  near  future  we  hope  to  report  an 
analysis  of  the  important  convolution  integral  approach. 
All  the  new  schemes  are  useful  in  studies  concerning  pos¬ 
sible  medical  effects  of  human  exposure  to  pulsed  electro¬ 
magnetic  fields.  This  type  of  application  requires  the  accu¬ 
rate  representation  of  amplitude  and  phase  information  for 
each  frequency  component  in  the  problem,  and  commonly 
involves  the  execution  of  thousands  of  timesteps.  Thus,  the 
numerical  method  of  choice  should  correctly  model  any  en¬ 
ergy  loss  due  to  physical  mechanisms,  and  it  should  intro¬ 
duce  as  little  phase  error  as  possible,  in  addition  to  being 


stable  for  all  possible  waves  that  can  be  supported  on  the 
computational  grid. 

The  stability  and  phcise  error  properties  of  the  standard 
FD-TD  approach  are  well  known  [2].  In  sharp  contrast, 
little  is  known  of  the  numerical  error  properties  of  the  dis¬ 
cretizations  of  the  extensions  of  FD-TD  to  dispersive  media 
modeled  with  ODEs  [3]-[5].  Questions  that  should  be  an¬ 
swered  before  such  extensions  are  used  in  production  codes 
include  the  following: 

1)  Is  the  stability  restriction  of  the  FD~TD  that  is  used  to 
model  the  PDE  part  of  the  problem  altered  by  appending  the 
discretized  (to  second  order  of  accuracy)  ODEs  that  model 
the  maternal  dispersion  ? 

2)  What  are  the  artificial  dissipative  properties  of  the  com¬ 
posite  scheme  ? 

3)  How  does  the  phase  error  now  depend  on  the  discretiza¬ 
tion  parameters  ? 

4)  Are  the  medium  constants  (relaxation  time,  resonance 
frequency,  etc.)  altered  by  the  discretization  ?  If  so,  how  ? 

We  answer  these  questions  by  performing  a  stability  and 
phase  error  analysis  [8]  of  the  coupled  PDE-ODE  difference 
system  that  models  one-dimensional  pulse  propagation  in 
linear  dispersive  media  of  Debye  [9]  and  Lorentz  [10]  type. 
The  analysis  is  novel,  and  its  purpose  is  to  provide  insight 
to  the  spurious  numerical  attributes  of  FD-TD  extensions 
for  dispersive  dielectrics.  We  give  many  details  for  the  one¬ 
dimensional  case  to  bring  forward  the  salient  features.  In 
Section  IV  we  show,  without  detailed  derivations  and  Fig¬ 
ures,  the  stability  polynomial  and  the  dispersion  relation 
of  the  extensions  studied  herein  for  two/three-dimensional 
propagation  in  Debye  and  Lorentz  dielectrics.  In  the  Dis¬ 
cussion  section  some  guidelines  will  be  provided  for  us¬ 
ing  the  schemes  given  the  medium  parameters.  The  main 
analysis  shows  that  for  accuracy  the  discretization  has  to 
adequately  sample  the  shortest  timescale  in  the  problem 
regardless  of  whether  it  is  the  incident  pulse  timescale, 
the  medium  relaxation  timescale,  or  the  medium  resonance 
timescale.  The  implication  of  this  is  that  one  usually  will 
have  to  resolve  a  scale  that  may  be  a  tenth  (usually  a  hun¬ 
dredth  or  a  thousandth)  of  that  of  the  incident  pulse.  Such 
disparity  of  scales  occurs  often  in  practice  where  experi¬ 
mental  data  indicates  that  typical  medium  timescales  are 
of  0(10“®  —  10“^^)  seconds  while  the  duration  of  pulses 
of  large  amplitude  (of  interest  to  health  and  safety  ana¬ 
lysts)  is  of  the  order  of  hundreds  of  nanoseconds.  Also, 


the  Courant  number  should  be  chosen  close  to  its  max¬ 
imum  value  for  stability  so  that  the  phase  error  is  min¬ 
imum  for  a  given  timestep.  Another  important  finding 
is  that  the  discretization  in  [3]  for  a  Lorentz  medium  is 
weakly  unstable  for  wavenumbers  k  such  that  AAz  >  7r/2. 
These  wavenumbers  correspond  to  unresolved  waves  (less 
than  4  cells  per  wavelength)  and  numerical  noise  (grid-scale 
oscillations).  Due  to  the  instability  such  low-amplitude 
waves  may  contaminate  well-resolved  events  in  computa¬ 
tions  encompassing  thousands  of  timesteps,  a  case  typical 
of  computations  in  dispersive  media  which  usually  involve 
0(10^-10*)  timesteps.  The  discovered  instability  will  tend 
to  slowly  amplify  noise  instead  of  producing  detectable  nu¬ 
merical  overflow.  The  Lorentz  scheme  in  [5]  is  found  to  be 
stable  for  0  <  ifeAx  <  tt,  and  to  slightly  damp  the  range 
JfeAx  >  7r/2  so  that  the  timestepping  itself  will  remove  the 
numerical  noise  instead  of  amplifying  it. 


is  determined  by  |^|.  Substituting  (2)  in  (1)  we  collect  the 
coefficients  of  the  vector  x  in  the  form  of  a  homogeneous 
linear  system  Ax  =  0  where  the  matrix  A  involves  the 
discretization,  and  the  medium  parameters.  Then  we  seek 
non-zero  solutions  for  x  by  setting  (iei(A)  =  0.  Extensive 
manipulations  of  the  entries  in  A  result  in  an  algebraic 
polynomial  equation  for  The  polynomial,  whose  solu¬ 
tions  give  ^  as  a  function  of  the  medium  parameters,  the 
timestep,  and  the  quantity  k/^  (—.  27r/^/pjyi0,  ——points 
per  wavelength),  is 


3  p^(/t  +  2)-6foo-/tC,  2 

^  ^  2eoo+h€.  ^ 

p\h-2)  +  5€oo-h€, 

2coo  4"  k,€g 

26oo  -  h€g  _  Q 

2^00  4" 


(3) 


II.  Stability  Analysis 
A,  Debye  Medium 

For  brevity  we  will  not  include  details  of  derivations  since 
they  are  quite  lengthy.  Rather,  we  will  describe  any  miss¬ 
ing  steps  so  the  reader  can  replicate  them.  Throughout,  we 
assume  the  reader  will  refer  to  the  appropriate  reference  for 
the  relevant  physical  equations.  For  the  sake  of  complete¬ 
ness  we  include  the  difference  schemes  for  the  approaches 
treated  herein. 

We  begin  by  considering  the  scheme  presented  in  [3]  for 
Debye  dispersive  media.  The  difference  equations  for  the 
update  of  the  magnetic,  displacement,  and  electric  fields 
read: 


ri  }  f]  ’ 

2r€oo-€gAt  ^ 


where  n  is  the  discrete  time  index,  j  is  the  discrete  spatial 
index,  rj  =  2r€oo  4-  c.  At,  r  is  the  medium  relaxation  time, 
€oo  is  the  infinite  frequency  permittivity,  e,  is  the  static 
permittivity,  and  /lo  is  the  permeability  of  the  vacuum. 
Neglecting  the  effects  of  boundary  conditions  and  initial 
conditions  a  solution  of  (1)  is  defined  as  follows 


Hf  ^ 

DJ  = 


Cn^ikjA 


In  (2)  the  complex  valued  vector  x  =  d,  e}^  is  the 
eigenvector  of  the  difference  system,  (  is  the  complex  time- 
eigenvalue  we  wish  to  find  and  whose  magnitude  will  deter¬ 
mine  the  stability  and  dissipation  properties  of  the  differ¬ 
ence  equations,  and  k  is  the  real  wavenumber  of  the  arbi¬ 
trary  harmonic  wave  component  whose  stability  and  decay 


where  p  =  2i/sin^,  h  =  At/r,  v  =  and  all  per- 

mittivities  are  now  relative  to  The  speed  Cq©  Is  the 
maximum  wavespeed  in  the  problem  and  is  given  by  Cq©  = 
c/y^e^.  The  speed  of  light  in  free-space  is  c. 

Next  we  considered  the  scheme  presented  in  [4].  The 
analysis  determines  that  the  two  schemes  are  identical  in 
terms  of  their  stability  polynomial  and  phase  error  (see 
Section  ///.A),  provided  we  make  the  following  identifica¬ 
tions  and  definitions  for  the  medium  parameters  used  in  (3) 
above,  i.e.,  set  foo  =  1  and  e,  =  1+X  where  x  is  the  suscep¬ 
tibility  used  in  [4].  Taking  r  =  8.1  x  10"^^sec,  e,  =  78.2, 
and  foo  =  li  we  calculate  numerically  with  IMSL  routines 
the  3  roots  of  (3)  for  a  variety  of  h  and  Courant  numbers 
u  as  functions  of  JkA.  These  medium  parameters  describe 
the  main  relaxation  of  water  in  the  microwave  range  of 
frequencies. 

By  examining  the  root  of  largest  magnitude,  77uix|^|,  as 
a  function  of  JfeA  in  the  range  0  <  iA  <  tt  we  determined 
that  the  difference  equations  are  stable  since  it  was  always 
that  max\C\  <  1  for  any  kA  in  the  range  considered  as  long 
as  1/  <  1  and  h  arbitrary.  In  such  a  graph  the  amount  of 
artificial  dissipation  introduced  by  the  differenced  ODEs 
is  evident  since  it  is  known  that  the  FD-TD  differenced 
PDEs  should  always  have  |  =  1  for  all  A: A  when  v  <  1. 
This  is  also  obtained  from  (3)  by  solving  for  ^  with  /i  =  0. 
The  result  is  shown  in  Fig.  1  for  i/  =  1.  In  the  Figure 
we  consider  the  values  h  =  0.1,  =  0.01,  and  h  =  0.001 

corresponding  respectively  to  a  resolution  of  10,  100,  and 
1000  timesteps  per  relaxation  time,  t.  It  turned  out  that 
Tnax\^\  >  1  (instability)  whenever  i/  >  1  (or  i/  > 
if  Coo  <  c)  regardless  of  the  medium  parameters.  Thus, 
the  one-dimensional  stability  restriction  of  the  standard 
FD-TD  scheme  is  preserved  by  both  extensions  for  Debye 
media.  Further,  since  at  a  given  space-time  point  the  dis¬ 
persion  is  the  material’s  memory  of  past  Electric  field  val¬ 
ues,  found  along  the  axis  of  the  light-cone  extending  into 
the  past  with  vertex  at  the  given  spatial  point,  the  sta¬ 
bility  restriction  of  the  standard  FD-TD  in  two  and  three 
dimensions  will  also  be  preserved  by  these  schemes.  This 
is  because  the  discrete  dispersion  variables  are  appropri- 
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ately  centered  with  respect  to  the  discrete  Electric  field 
so  for  0  <  v  <  l/VSim,  where  dim(=  1,  2,  3)  is  the  spa¬ 
tial  dimension,  the  characteristics  of  the  overall  scheme  are 
the  same  as  those  of  the  standard  FD-TD  and  include  the 
physical  characteristics  whose  slope  is  icoo- 

B.  Lorentz  Medium 

For  the  Lorents  medium  case,  the  schemes  of  [3]  (here 
called  the  JHT  scheme),  and  [5]  (here  called  the  KF  scheme), 
differ,  both  in  the  structure  of  their  difference  equations, 
and  in  their  stability  and  phase  error  properties. 

We  begin  with  the  scheme  given  in  [3]  for  Lorentz  media. 

The  difference  equations  now  are  as  in  (1),  above,  for 

and  but  the  electric  field  update  is  performed  with 

“  3  ’  # 

=  -E^  +  -E2~^  +  +  --D”  + 

J  rj  7]  7] 

where  now  r,  =  wjAi’e.  -1-  25At£oo  +  2£oo,  the  other  coefB- 
cients  being 


P  =  -(w’ At’cj  -  2JAtCoo  +  2coo) 

7  =  w^At^  +  25  At  +  2  (5) 

C  =  -4. 

0  =  wlAt^  -25At  +  2. 

In  (5)  Uo  is  the  medium  resonance  frequency,  5  =  l/2r  with 
T  being  the  relaxation  time,  and  the  rest  of  the  parameters 
are  as  above  for  the  Debye  medium.  Substituting  (2)  into 
the  difference  equations  we  obtmn  the  associated  stability 
polynomial,  now  of  fourth  order  in  (  because  (4)  involves 
an  extra  level  of  storage, 

+  L  =  0  (6) 

n  T1 


to  the  medium  resonant  frequency  so  is  a  measure  of  how 
well  the  resonant  period  of  the  medium  is  resolved  by  the 
time  discretization. 

Next,  we  analyze  the  scheme  of  [5].  The  difference  equa- 
tlons  given  there  are 


+  P"/«o 


7  =  2  +  hi  + 


fl  =-4 


xB-*‘ + T'Vr  = 

where  x  =  l"  =  1/’’.  free-space  permittivity, 

and  w,  is  the  plasma  frequency.  In  terms  of  the  pwame- 
ter  definitions  used  for  the  Lorentz  scheme  of  [3)  it  is  that 
_  1).  This  scheme  is  nice  for  the  following  rea¬ 
son:  The  second,  third,  and  fourth  difference  equations  in 
(8)  arc  implicit  in  the  variables  P"'*’*!  while 

the  first  and  second  are  explicit  in  As  a  re¬ 

sult,  the  stability  properties  are  independent  of  the  iiiedium 
time  constants,  and  the  Courant  restriction  for  stability  of 
the  standard  scheme  will  be  retained  even  in  two  and  three 
dimensions.  The  partial  implicitness  in  (8)  is  not  trou¬ 
blesome  since  the  matrix  relating  Py  ^ 
inverted  just  once,  except  when  At,x,  wj,  and  r  are  func¬ 
tions  of  time.  A  solution  of  (8)  is  of  the  form 


Substituting  (9)  in  (8)  and  carrying  out  the  steps  indicated 
above  we  arrive  at  the  stability  polynomial  for  this  scheme 


for  which  now 


^  =  2  --  hi  -|-  47r^h3 

a  =  — 8foo  ~  2/ii€oo  ~  87r^h3£<  C^) 

/3’  =z  12«oo  +  iv^hle. 

5  —  "Scoo  “h  2hi€oo 

rj'  =  2«oo  -  hi«oo  +  4e,x*h3, 

where  now  hi  =  At/r,  hi  =  At/To  (T<,  =  2ir/u)o),  and  t]  — 
2£oo  +  hieoo  +4£.7r’h’.  The  period  To  is  that  corresponding 


7  =  4  -h  2hi  +  Air^h\ 

0  +  iir^hl 

C'  =  4  -  2hi  +  Ait^h] 

a  =  -16  -  4hi 

/3'  =  24  -  8<,7r^h| 

s'  =  -16  +  4hi 

ij  =  4  —  2hi  +  4£,n’^h2, 


and  7]  =  A-\‘2hi-\- Ae^TT^ hi-  Equation  (10)  is  of  4th  order  in 
(  not  because  an  additional  time  level  was  introduced  but 
because  the  second  order  in  time  ODE  for  the  polari2ation 
was  transformed  to  a  first  order  system  in  accordance  with 
Maxwell  equations. 

The  parameters  used  in  the  two  approaches  considered 
will  address  the  same  physical  problem  if  we  set  x  = 
and  Wp  =  ”  1)  parameters  used  in  [5]  using 

values  for  Wo,  r,  Cooi  and  e,  from  the  numerical  example  of 
[3],  With  these  values  we  study  the  4  roots  of  the  stability 
polynomials  (6)  and  (10)  using  the  parameter  definitions 
(8)  and  (11)  for  each  respectively.  The  values  used  here 
are  r  =  1.786  x  10"^®sec,  Coo  =  1,  e.  =  2.25,  and  Uo  = 
4xl0^®rad/sec.  The  results  are  shown  in  Fig.  3a)  for  i/  =  1 
and  hi  =  0.1,  and  in  Fig.  3b)  for  i/  =  1  and  hi  =  0.01. 
For  the  scheme  of  [5]  it  is  determined  that  max\^\  >  1, 
which  corresponds  to  instability,  occurs  only  when  i/  >  1 
regardless  of  the  medium  parameters.  However,  for  the 
scheme  ^  [3]  we  observe  instability  for  the  range  AA  >  Tt/2 
and  some  ratios  hi  =  At/r. 


III.  Phase  Error  Analysis 
Now  we  will  determine  the  phase  error  of  each  approach 
considered  in  the  previous  section.  After  proper  parameter 
identifications  the  modeling  approaches  of  [3]  and  [4]  for 
Debye  media  have  identical  phase  error  properties.  We 
will  consider  the  following  definition  of  phase  error 


l.D,L 

Kex 


(a;) 


(12) 


for  real  w,  where  is  the  dispersion  relation  of  the 

continuous  models  of  Debye  (D)  and  Lorentz  (L)  media 
given  by 


+  2iS(j 
u)^  —  2iSD  ' 


(13) 


and  (orAt)  will  be  determined  in  the  subsections  be¬ 
low.  Note  that  for  both  types  of  media  we  have  set  Coo  = 
so  the  expressions  simplify  (also  this  value  is  used  in  the 
numerical  experiments  presented  in  the  references).  We 
have  used  =  ^oXt  X  —  ”  1/r  to  relate  the 

Lorentz  mediums  of  [3]  and  [5]‘.  The  uf  range  considered  is 
such  that  0  <  wAt  <  ir  given  At.  The  spatial  step  A  is 
determined  from  At  and  i/. 


A.  Debye  Medium 

Only  the  scheme  of  [3]  will  be  considered  since  it  is  iden¬ 
tical  to  that  of  [4]  for  Coo  =  To  determine  the  numerical 
dispersion  relation  we  now  substitute 


(14) 


in  the  difference  equations  (1),  with  (h,d,  e)^  being  an  ar¬ 
bitrary  complex  vector,  and  manipulate  the  result  so  that 
again  we  obtain  a  homogeneous  linear  system  whose  de¬ 
terminant  of  the  coefficient  matrix  we  demand  to  vanish. 
The  dependence  of  on  the  frequency  cj  and  on  the 

various  medium  and  discretization  parameters  is  thus  ob¬ 
tained.  The  result  is 


lD 


.(wAt)  =  -  sin 


cos 


— 


1  wAt 

-cos  — 


—  VjJSu 


1.  (15) 


where  s„  =  ■  By  inspecting  (15),  and  comparing 

it  to  given  above,  two  features  emerge  that  are  solely 
due  to  the  discretization  of  the  ODE  involved.  The  re¬ 
laxation  time  r  of  the  medium  is  now  Tnum  = 
i.e.,  the  medium  actually  modeled  by  the  numerics  is  one 
with  higher  relaxation  time  constant.  This  is  the  source  of 
the  artificial  dissipation  exhibited  by  the  maximum  root, 
max|^l,  of  (3).  Further,  since  cos(a;At/2)  0  for  a; 
7r/At,  we  infer  that  the  corresponding  frequency  compo¬ 
nent  will  be  adversely  affected  by  the  scheme  since  it  will 
experience  a  medium  of  unphysically  high  constant  conduc¬ 
tivity  equal  to  Co(ft  -  l)/rnum-  Such  artificial  dissipation 
can  be  controlled  by  choosing  Ai  so  that  cos  cj At/2  ^  1 
across  the  range  of  frequencies  present  in  the  short-pulse 
that  propagates  in  the  medium.  For  i/  =  1,  Fig.  2  shows 
the  dependence  of  the  phase  error  (12)  on  the  number  of 
timesteps  per  relaxation  time,  h  =  At/r,  ets  a  function  of 
luAt. 


B,  Lorentz  Medium 

For  the  scheme  in  [3]  pertaining  to  Lorentz  media,  we 
follow  the  same  procedure  as  with  the  Debye  medium  but 
now  the  third  difference  equation  in  (1)  is  replaced  with  (4). 
Using  (14)  we  obtain  for  the  numerical  dispersion  relation 


k^ 


r  .  .  2  .  ifWA 

sm" 


cos  uAt  +  iSu 
u'^sl  -w^coswAf + 


Equation  (16)  shows  how  this  extension  of  FD-TD  misrep¬ 
resents  the  relaxation  time  by  enforcing  < 

and  the  resonance  frequency  of  the  medium  by  enforcing 
^num  _  a;o>/cos(cjAt). 

Next,  we  determine  the  phase  error  of  the  scheme  in  [5]. 
Substituting 


(17) 


in  the  difference  equations  (10)  we  obtain  the  dispersion 
relation 


^num(‘^At)  = 


2  .  _i,w  A 
—  sin  — —3u 
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ju^sl  -  e,ul  cos^  ^  +  2iScoa 
Y  cos^  +  2i6  cos 

now  for  7  of  [5]  equal  to  2S  (S  is  used  in  [3]).  From  (18) 
we  see  that  this  approach  also  misrepresents  the  relax¬ 
ation  time  by  enforcing 

same  as  the  misrepresentation  in  the  Debye  medium  case. 
The  numerical  resonance  frequency  of  the  medium  is  now 
^»«m  _  £,;^cOs(wAi/2). 

^Phase  error  comparisons  between  the  two  approaches  us¬ 
ing  (12)  are  shown  in  Figs.  4a)  and  4b)  for  u  =  I,  and  ratios 
hi  =  At/T=  0.1  and  0.01  respectively. 


IV.  2/3  Dimensional  Extensions 
Herein  we  assume  the  reader  is  familiar  with  the  two  dimen¬ 
sional  Transverse  Magnetic  polarization  formulation  of  the 
standard  FD-TD  [1].  The  unknown  two-dimensional  fields 
are  H,,  The  discretized  dispersion  ODE  of 

all  schemes  is  unchanged  in  higher  dimensions.  We  set 
A®  =  Ay  =  A. 

To  determine  the  stability  for  the  Debye  extension  of  the 
scheme  of  [3]  we  proceed  as  in  Section  / 1  .A.  The  difference 
equations  now  reflect  the  existence  of  the  extra  fields  and 
the  two  spatial  variables.  Accordingly,  (2)  is  extended  as 
follows 


f  H?(i.m) 
1 

]  D:u,m) 

[  K{j,m) 


f  1 

^  ^  I  Jk  coi^+mfc  iin^)  A  ^  (19) 

i- 

I  C,  j 


where  (j,  m)  are  the  discrete  spatial  indices,  k  —  yjk^  -f  kj 
is  the  wavenumber  magnitude,  and  9  is  the  angle  of  prop¬ 
agation  with  respect  to  the  grid  ®— axis.  Equation  (3)  be¬ 
comes 

p^(2  -I-  h)  -  8foo  -  2he.  3  12foo  -4p^  3 
^  ^  2eoo  +  he,  2eoo  +  he. 


+ 

•+ 


p^(2  -  h)  -  8eoo  +  2he. 

2eoo  +  he, 

2eoo  —  he. 


2eoo  -k  he. 


=  0, 


(20) 


in  the  difference  equations  and  proceeding  as  in  one  di¬ 
mension.  After  performing  the  calculations  we  find  that 
the  dispersion  relation  for  real  frequency  u  is 

sin^  p(w,  9,  A)  -f-  sin^  C(w,  9,  A)  = 


2  Vcosiiiyt  -ie„USu 
ic03!=^-iw5„ 

r  2 


(22) 


where  p  =  k^^,j,{o> £^t)  cos  0( A/2),  ^  sin  9[A/2) 

Again,  for  too  =  1  the  two-dimensional  extension  of  the  De¬ 
bye  scheme  of  [4]  is  again  identical  to  that  of  [3). 

The  stability  of  the  Lorentz  scheme  of  [5]  can  be  ad¬ 
dressed  with  the  discussions  in  the  end  of  Section  1 1.  A 
and  below  equation  (8).  The  relevant  polynomial  for  the 
Lorentz  scheme  of  [3]  is 

^  n  r, 

p2(0-C)-l-3a  +  3/3  2 
rt 

a-3/3-_gp^  (23) 

V  ^ 

and  the  various  symbol  definitions  can  all  be  found  in  equa¬ 
tions  (4)  and  (5).  In  (23)  p^  is  the  same  as  for  equation 
(20).  The  numerical  dispersion  relations  for  all  the  Lorentz 
medium  schemes  are  easily  obtained  from  equations  (16) 
and  (18)  by  inspection  considering  that  dispersion  is  a  tem¬ 
poral  phenomenon  and  does  not  enter  in  geometrical  con¬ 
siderations.  As  a  result,  the  two-dimensional  dispersion 
relations  are 


sin*  p{u,  9,  A)  -t-  sin*  C(‘*'i  A)  —  A*,  (24) 

where  A  is  the  argument  of  the  sin“^  function  in  (16)  and 
(18),  and  p  and  <  are  defined  as  for  (22)  but  with  ^num 
replacing  Notice  that  (22)  can  also  be  derived  m 

the  same  way  with  A  now  being  the  argument  of  the  sin 
function  in  (15). 

In  three  dimensions  the  dispersion  relation  for  all  schemes 
(Debye  and  Lorentz)  considered  herein  will  be 

sin*  p(w,  9,  A)  +  sin*  C(w,  9,  V>,  A) 

-J-  sin*  7r(a;,  0, 1^,  A)  =  A*.  (25) 


where  now  p*  =  4i/*(sin*  p  +  sin*  ^),  p  —  kcos6(A/2), 
and  (  =  ksin9(A/2).  The  numerical  dispersion  relation  is 
obtained  by  using 


'  H»(i,m) 

I  E2(j,m) 


> 


1 

^  I  gi(OJk?„„eoia+mif„„finfl)A-unA(] 


(21) 


In  (25),  A  is  as  above,  p  =  k^^^{it>At)  cos  9  sin  V’(A/2),  C  — 
Jkf,;m(wAt)sin0sin^(A/2),  tt  =  k°i^iu)At)co3rl,{A/2), 
with  ^  being  the  angle  of  propagation  of  the  harmonic  wave 
with  respect  to  the  z-axis  and  now  9  being  the  angle  be¬ 
tween  the  grid  ®-axis  and  the  projection  of  the  wavevector, 
k^^^{wAt)  =  xk^'^  +  yk^'’^  +  0“  the  ®  -  y  plane. 

V.  Discussion 

The  finite  difference  approaches  of  [3]  and  [4]  for  Debye  me¬ 
dia  have  been  determined  to  be  equivalent  in  any  number 
of  dimensions.  Therefore,  we  studied  the  stability  polyno¬ 
mial  (3)  and  the  phase  error  (15),  obtained  for  the  scheme 


I  e. 


of  [3]  in  one  dimension.  Fig.  1  shows  the  magnitude  of 
the  largest  of  the  3  roots  of  the  stability  polynomial  as  a 
function  of  the  product  which  may  be  thought  of  as 
wavenumber  if  A  is  fixed,  or  as  inverse  of  points  per  wave¬ 
length  (equivalently,  as  A)  if  k  is  viewed  as  fixed.  While 
the  Figure  shows  all  roots  to  be  inside  the  unit  circle  (hence 
the  scheme  to  be  stable),  it  also  indicates  that  the  differ¬ 
ence  equations  arc  now  spuriously  dissipative.  This  is  in 
contrast  to  the  case  where  FD-TD  solves  Maxwell’s  equa¬ 
tions  in  free-space  (the  solid  line  in  these  Figures),  or  when 
the  dispersion  ODEs  are  solved  separately  from  Maxwell’s 
equations,  i.e.,  set  =  0  in  (3)  and  solve  for  the  roots.  It 
turned  out  that  for  i/  <  1  it  was  max\i\  <  1  always,  and 
that  T7wx|^|  >  1  only  for  i/  >  1.  Therefore  these  exten¬ 
sions  of  FD-TD  for  Debye  media  preserve  the  stability  re¬ 
striction  of  the  standard  FD-TD  scheme  in  non-dispersive 
dielectrics.  The  graph  also  shows  that  the  artificial  dissi¬ 
pation  is  strongly  dependent  on  the  quantity  A.  We  see 
that  h  has  to  be  sufficiently  small  so  that  not  much  en¬ 
ergy  is  lost  after  the  large  amount  of  timesteps  usually 
needed  to  propagate  short  pulses  any  appreciable  distance 
into  a  medium.  One  guideline  for  a  user  would  be  that  the 
timestep  chosen  must  be  at  least  of  O(10~^r),  where  r  is 
the  relaxation  time  in  the  Debye  medium,  or  the  smallest 
relaxation  time  in  a  medium  modeled  with  multiple  Debye 
relaxation  mechanisms.  In  Fig.  2  the  phase  error  (12)  is 
graphed  as  a  function  of  uAt  for  the  same  three  values  of 

h.  The  horizontal  axis  again  has  a  triple  interpretation, 

i. e,  frequency,  timestep,  or  points  per  period.  Again,  we 
see  that  At  --  O(10'“^r)  is  optimal.  In  both  Figures  i/  =  1, 
i.e.,  the  maximum  Courant  number  for  stability.  It  turned 
out  that  if  i/  <  1  then,  for  a  fixed  accuracy,  the  timestep  re¬ 
striction  with  regard  to  the  relaxation  time  is  more  severe, 
thus  another  guideline  is  that  these  schemes  should  be  run 
close  to  their  Courant  stability  limit,  and  this  will  result  in 
minimal  run  time  since  the  timestep  will  then  be  optimally 
maximum.  One  must  point  out  that  these  results  do  not 
depend  on  the  points  per  wavelength  used  to  sample  the 
incident/propagating  pulse.  If  the  pulse  duration  is  longer 
than  the  relaxation  time,  one  must  still  maintain  a  sam¬ 
pling  rate  that  resolves  the  shortest  timescale,  which  is  the 
relaxation  time  in  this  case.  If  the  incident  pulse  duration 
is  shorter  than  r  one  decides  on  the  sampling  in  the  usual 
way  as  long  as  =  10“^. 

Figs.  3a)  and  3b)  are  concerned  with  stability  compar¬ 
isons  of  the  JHT  scheme  and  the  KF  scheme,  for  Lorentz 
media.  We  see  that  the  JHT  scheme  is  unstable  (mai|^|  > 
1)  for  kA  >  7r/2  at  i/  =  1  if  /ii  is  not  chosen  correctly, 
whereas  the  standard  FD-TD  and  the  KF  scheme  are  sta¬ 
ble  for  all  wavenumbers.  It  turns  out  that  the  JHT  ex¬ 
tension  has  a  stability  limit  that  is  more  restrictive  than 
that  of  FD-TD  in  non-dispersive  dielectrics.  Reducing  hi 
to  at  least  10“^  stabilizes  the  scheme.  We  attribute  this 
instability  to  incorrect  time  centering  of  the  local-in-time 
constitutive  relation  with  respect  to  the  time  centering  of 
the  discretized  PDEs.  From  the  difference  equation  (4)  for 
the  constitutive  relation  it  is  seen  that  dD/dt^  d^D/di?  are 
centered  at  the  Tvth  timestep  while  in  the  Maxwell’s  equa¬ 


tions  discretization  dDfdt  is  centered  at  the  timestep. 
This  is  also  indicated  by  the  presence  of  sinusoidal  factors 
with  argument  a;At  and  ufAt/2  in  the  dispersion  relation 
for  the  scheme  (16).  On  the  other  hand,  the  KF  scheme 
centers  the  ODE  part  correctly  and  does  not  exhibit  this 
instability.  The  phase  error  of  the  two  Lorentz  schemes  are 
shown  in  Figs.  4a)  and  4b)  for  i/  =  1  and  two  values  of 
hi.  Running  the  schemes  with  i/  <  1  is  again  discouraged. 
For  these  schemes  we  notice  that  it  is  enough  to  choose 
At  ~  O(10”^r),  where  r  =  1/(2J).  Note  that  the  KF 
scheme  exhibits  consistently  less  phase  error  than  the  JHT 
scheme.  We  have  not  investigated  the  effect  of  varying  hj 
since  for  the  numerical  values  of  the  medium  parameters 
used  here,  r  ^  To . 

In  Section  IV  we  presented  the  stability  polynomials  of 
the  Debye  medium  schemes  and  of  the  weakly  unstable 
Lorentz  scheme  of  [3]  for  two-dimensional  TM  propagation. 
All  schemes  that  correctly  center  the  dispersion  variables 
preserve  the  known  stability  restriction  of  the  standard  FD- 
TD  in  higher  dimensions.  Numerically  solving  for  the  roots 
of  these  polynomials  and  determining  the  largest  root,  as 
a  function  of  the  discretization , parameters  for  each  angle 
of  propagation  of  the  harmonic  wave  given  the  medium 
parameters,  one  can  easily  ascertain  the  degree  of  decay 
each  mode  will  sustain  and  thus  will  be  guided  towards  the 
selection  of  optimum  parameters.  The  Debye  schemes  and 
the  Lorentz  scheme  of  [5]  are  stable  for  u  <  l/\/2  regardless 
of  the  value  of  /i,  or  hi.  The  Lorentz  scheme  of  [3]  is  again 
weakly  unstable.  The  dispersion  relations  of  all  schemes 
are  easy  to  deduce  in  higher  dimensions,  once  their  one¬ 
dimensional  form  is  known.  We  give  all  dispersion  relations 
in  two  and  three  dimensions.  Now  the  wavenumber  cannot 
be  explicitly  determined  as  a  function  of  the  real  frequency 
u).  An  iterative  scheme  must  be  used  to  calculate  k^^^ 
for  each  w,  0,  given  the  medium  parameters,  using  the 
non-dispersive  value  as  a  starting  value. 

Finally,  we  simulated  the  propagation  of  a  1  nanosec¬ 
ond  duration  square-modulated  sine  wave  with  carrier  fre¬ 
quency  at  10  GHz  normally  incident  on  a  Debye  medium 
(fj  =  80.35,  eoo  =  LO,  r  =  8.13  picoseconds)  half-space 
from  the  air  side.  Six  runs  were  made  with  a  code  written 
from  reference  [3],  The  timestep  was  successively  halved 
(starting  with  h  ^  0.05)  by  halving  h  for  fixed  r  and 

=  1.  A  time  trace  of  the  electric  field  at  a  depth  of  10  mm 
was  recorded  with  a  window  of  2  nanoseconds.  The  num¬ 
ber  of  timesteps  executed  in  each  run  was  2  x  10~®/ At 
0(10^  —  10^).  Fig.  5  shows  the  time  traces  from  the  first 
three  runs  and  how  they  converge  with  the  reduction  of 
the  timestep.  Further  reduction  shows  that  the  timestep 
for  which  h  =  0.00615  is  adequate  for  a  converged  re¬ 
sult.  In  Fig.  6,  the  h  =  0.00615  result  is  indistinguishable 
from  the  h  —  0.003075  result.  Graphs  of  the  dissipation 
and  phase  errors  (not  shown  here)  reflect  the  fact  that, 
in  terms  of  numerical  errors,  the  computed  results  should 
not  change  much  with  further  timestep  reduction  past  the 
value  for  which  h  =  0.00615.  With  only  this  information  we 
could  have  avoided  the  six  convergence  runs.  The  guideline 
At  --  (^(IO-^t)  is  confirmed. 


VI.  Summary 

In  this  paper  we  have  analyzed  four  FD-TD  based  schemes  ,c«ttering  in 

for  modeling  pulse  propagation  in  dispersive  niedia  of  De-  .  member  of  Ta«  Bd. 
bye  and  Lorentz  type.  The  medium  dispersion  is  described 
by  ODEs  representing  the  local-in-time  constitutive  rela- 
tion  f3]  or  the  dynamic  equation  for  the  polarization  evo¬ 
lution  driven  by  the  electric  field  [4]-[5].  We  found  that 
the  extensions  considered  here  do  not  preserve  the  non- 
dissipativc  character  of  the  standard  FD-TD  which  is  good 
for  tracking  individual  frequency  components  for  a  large 
amount  of  timesteps,  and  that  the  ease  of  limiting  artifi¬ 
cial  dispersion  is  lost  since  now  there  are  more  parameters, 
in  addition  to  i/,  that  control  it.  The  overall  difference  sys¬ 
tems  are  more  dispersive  than  the  standard  FD-TD,  and 
their  accuracy  depends  strongly  on  b#>w  well  the  chosen 
timestep  resolves  the  shortest  timescales.  All  parameters 
which  play  a  major  role  in  the  elimination  of  artificial  dis¬ 
sipation  and  phase  error  for  these  schemes  have  been  iden¬ 
tified.  The  higher  dimensional  form  of  our  results  has  also 
been  given. 
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1.  The  maximum  stability  eigenvalue  of  the  Debye  medium  difference  schemes  versus 
A:A(=  iTzfNppJ)  for  u  =  CgoAt/A  =  1.0. 
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FIGURE 


3.  Comparison  of  the  stability  eigenvalues  of  the  JHT  and  KF  Lorentz  medium  differ¬ 
ence  schemes  for  a)  hi  =  0.1,  and  b)  hi  =  0.01. 
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3.  CompaiUon  of  the  stebiUty  eigenvelees  of  the  JHT  and  KF  Lorenta  medium  differ- 
ence  schemes  for  a)  hi  =  0.1,  and  b)  hi  =  0.01. 
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4.  Compaxison  of  the  phase  errors  attained  by  the  JHT  and  KF  Lorentz  medium 
difference  schemes  for  a)  hi  =  0.1,  and  b)  hi  =  0.01.  Note  the  difference  in  scales 
between  the  two  Figures. 


4.  Comparison  of  the  phase  errors  attained  by  the  JHT  and  KF  Lorentz  medium 
difference  schemes  for  a)  hj  =  0.1,  and  b)  hi  =  0.01.  Note  the  difference  in  scales 
between  the  two  Figures. 
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5.  Electric  field  computed  with.  FD-TD  for  Debye  medium  as  a  function  of  A  =  At/r. 
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6.  Same  as  Fig.  5  for  further  reduction  of  h.  Convergence  is  achieved  at  h  =  0.00615. 
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Th6  ^Vave  Hiararchy  for  Propagation  in  Rolaxing  Dioloctrics 


« 


Abstract 

We  consider  the  propagation  of  arbitrary  electromagnetic  pulses  in  anomalously 
dispersive  dielectrics  characterized  by  M  relaxation  processes.  A  partial  differential 
equation  for  the  electric  field  in  the  dielectric  is  derived  and  analyzed.  This  single 
equation  describes  a  hierarchy  of  M  +  1  wave  types,  each  type  characterized  by  an  at- 
tenuation  coefficient  and  a  wave  speed.  Our  analysis  identifies  a  “skin-depth  where  the 
pulse  response  is  described  by  a  telegrapher’s  equation  with  smoothing  terms,  travels 
with  the  wavefront  speed,  and  decays  exponentially.  Past  this  shallow  depth  we  show 
that  the  pulse  response  is  described  by  a  weakly  dispersive  advection-diffusion  equa¬ 
tion,  travels  with  the  sub-characteristic  advection  speed  equal  to  the  zero-frequency 
phase  velocity  in  the  dielectric,  and  decays  algebraically.  The  analysis  is  verified  with 
a  numerical  simulation.  The  relevance  of  our  results  to  the  development  of  numerical 
methods  for  such  problems  is  discussed. 
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1.  Introduction 


The  application  of  ultra-short  pulsed  fields  in  the  areas  of  radar,  hyperthermia,  and  bio¬ 
logical/environmental  imaging  is  imminent.  For  that  reason  there  is  a  need  for  a  thorough 
understanding  of  the  short-pulse  response  of  media,  such  as  biological  tissues,  soils,  humid 
atmospheres,  and  radar  absorbing  materials,  whose  dielectric  properties  are  described  by 
frequency-dependent  permittivity  models  fitted  to  experimental  data.  In  addition,  a  future 
extension  of  the  IEEE  C95. 1-1991  RF  exposure  standard  to  regulate  pulsed  fields  will  also 
require  qualitative  and  qiftintitative  understanding  of  this  sort  to  be  developed.  Since  the 
zdternative  to  actual  measurement  of  the  response  is  numerical  simulation  we  are  eJso  inter¬ 
ested  in  obtaining  anadytical  results  to  guide  the  development  of  robust  numericad  techniques 
for  such  applications  (see  Section  6). 

In  this  paper  we  consider  time-domaiin  electromagnetic  pulse  propagation  in  general  relzix- 
ing  dielectrics,  a  problem  relevant  to  studies  of  electromagnetic  interactions  in  homogeneous 
anomalously  dispersive  media.  We  will  be  concerned  primarily  with  biological  media  but  our 
results  can  be  used  directly  to  understand  pulse  propagation  in  other  dielectrics  provided 
their  permittivity  is  represented  by  relaocation  models.  The  Debye  model  [1]  (orientational 
relajcation  mechanism)  has  been  found  to  accurately  represent  the  experimentally  obtained 
frequency-dependent  permittivity  of  biological  media  [2]-[3],  of  homogeneous  dispersive  rock 
[4]-[5],  and  of  homogeneous  sea  ice  [6].  All  these  media  generadly  exhibit  several  relaxation 
mechanisms.  Accordingly,  experimentzdly  obtzdned  dielectric  data  is  fitted  to  the  complex 
relative  permittivity  function 


for  the  relevant  range  of  frequencies  uj  in  order  to  study  wave  propagation  in  a  given  mate¬ 
rial.  In  (1)  Coo  is  the  infinite-frequency  relative  permittivity,  e”  the  zero-frequency  relative 
permittivity  of  the  n-th  relaxation  mechanism,  and  the  n-th  relaxation  time  constant. 
The  phase  velocity  of  harmonic  waves  in  the  dielectric  is  =  cj ReyJ £(ca;),  with  c 

being  the  speed  of  light  in  vacuum.  Typical  pairameters  for  which  (1)  fits  water  data  in 
the  microwave  frequency  range  are  M  =  1,  e,  =  80.35,  too  =  1*0,  r  =  8.13  picoseconds 
[2].  Water  is  considered  here  since  it  is  a  major  constituent  of  human  tissue.  We  will  use 
the  inverse  Fourier  transform  with  (1)  in  the  electromagnetic  frequency-domain  constitutive 
relation,  D  =  tot[u))E  =  £0(^00^  +  Enii  ^l^se  the  system  of  Maxwell’s  equations  in 
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the  time-domain  with  M  ordinary  differential  equations,  one  for  each  relaxation  P'*  forced 
by  the  electric  field.  The  reformulation  in  Section  3  of  the  resulting  problem  (Equation  (3), 
Section  2)  as  a  signaling  problem  involving  a  single  partial  differential  equation  (p.d.e.)  for 
the  electric  field  E, 

=  (2) 

n=0 

facilitates  further  analysis.  The  coefficients  /3„  characterize  dissipative  effects  associated  with 
the  n— th  wave  type,  and  c„  is  the  corresponding  wave  speed.  Both  constants  are  related 
to  the  pzirameters  selected,  so  that  (1)  fits  the  experimental  data  for  a  given  medium.  All 
the  sub-characteristics  of  (2)  with  speeds  c„  (n  =  1, ...,  M)  he  between  the  characteristic 
with  speed  co,  and  the  zero  speed  characteristic  of  multiplicity  M.  Therefore,  the  relative 
permittivity  (1)  is  causal  and  its  rezd  and  imaginary  parts  are  related  through  the  Kramers- 
Kronig  relations  [7]. 

Equation  (2)  belongs  to  a  class  of  p.d.e.  which  describe  coexisting  waves  of  different 
types  [8],  e.g.,  non-dispersive,  dispersive,  diffusive.  In  Section  4  we  show  that  the  high- 
order  term,  n  =  0,  mzdnly  describes  the  propagation  of  the  penetrating  pulse  in  a  thin  layer 
near  the  <iir/di electric  interface.  In  this  layer  the  lower  order  terms  contribute  exponential 
decay  and  smoothing  of  the  pulse  envelope.  We  name  this  short  depth  the  skin-depth  for 
pulses.  The  pulse  propagates  there  with  speed  Cq  =  =  oo).  The  layer  s  width  is 

determined  by  the  time  constant  of  the  fastest  relaxation  since  it  will  be  the  first  one  to 
equilibrate,  i.e.,  ZMn  OCcor"”")  meters.  Each  lower  order  term,  1  <  n  <  M  -  1,  becomes 
sequentially  important  as  the  pulse  propagates  deeper  in  the  medium  while  the  remaining 
orders  introduce  smoothing,  dispersion,  and  diffusion.  Since  the  available  experimental  data 
shows  that  ^  1  typically  ,  the  low-order  term,  n  —  M,  will  describe  the  main  response 

which  will  travel  with  speed  cm  =  =  0)  in  the  bulk  of  the  medium  while  the 

higher  order  terms  will  introduce  high-order  dispersion  and  diffusion.  Thus,  the  response 
inside  the  dielectric  will  be  decaying  as  z"",  a  >  1/2.  The  achieved  value  of  a  depends  on 

the  pulse  frequency  content,  the  pulse  shape,  and  on  Af.  Explicit  solutions  will  be  given 

for  Af  =  1  whence  the  short-  and  long-time  response  respectively  satisfy  a  telegraphers 
equation  and  an  advection-diffusion  equation.  For  this  case  we  will  deduce  the  O(cot)  size 
of  the  “skin-depth”  for  pulses,  and  the  asymptotic  decay  rate  with  depth  of  the  peak 
electric  field,  or  with  time.  These  results  are  verified  with  a  numerical  simulation  in 

Section  5.  The  rate  of  decay  determined  herein  for  pulsed  electric  fields  may  be  of  interest  to 
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workers  in  the  fields  of  hyperthermia  or  ground  penetrating  radar.  The  signals  used  in  these 
technological  areas  are  continuous  waves  and  as  a  result  the  electric  field  in  the  dielectric 
(tissue  or  ground)  decays  exponentially  with  depth  [7]  due  to  the  imaginary  part  of  the 
frequency-dependent  permittivity.  Our  work  suggests  using  pulsed  electromagnetic  waves  to 
circumvent  shortcomings  aissociated  with  use  of  continuous  waves  in  such  applications. 


2.  Formulation 


A  Trzinsverse  Magnetic  pulse  is  obliquely  incident  on  a  homogeneous  dispersive  half-space 
from  the  air  side  (z  <  0)  at  an  angle  (f)ine  with  the  normal  to  the  air/dielectric  interface. 
The  dielectric  occupies  the  half-space  z  >  0.  The  equations  governing  the  scattering  and 
propagation  of  the  incident  pulse  Eire  the  time-dommn  Maxwell’s  equations  coupled,  thru 
the  electric  field,  to  M  ordinary  differential  equations  that  describe  the  evolution  of  M 
orientational  polarization  mechanisms  (P";  n  =  1, ...,  M,  see  (1))  of  Debye  type.  They  are: 

g/f,  dEy 

dt  ~  dz 


dH,  dEy 
dt  ~  dx 

dEy  dH,  dH, 

dt  ~  dz  dx  "“^1  dt 

^  -  P;);  n  =  1, M, 

ut  r” 


(3) 


where  60  and  fio  are  respectively  the  permittivity  and  permeability  of  the  vacuum,  and 
Ae”  =  —  600  •  Apart  from  and  [lo  3^1  other  parameters  in  (3)  are  obtained  by  fitting 

(1)  to  experimented  data  for  the  tissue  types  of  interest.  The  incident  electric  field  is  an 
arbitrary  plane  pulse  z,  t)  =  f{t  —  x  sin  ^incjc  —  z  cos  ^{nc/ c)  of  duration  Tp,  and  we 

assume  it  has  been  in  contact  with  the  interface  since  t  =  —  00.  Operational  considerations 
fix  the  incident  pulse  shape  /,  and  Tp.  On  the  interface,  z  =  0,  the  total  electric  field  is 
given  by  Ey{x,0,t)  =  f{t-x/v)  +  E‘^^^{x,0,t)  =  g{t-xlv),  where  t;  =  c/  sin  <f>inc,  and 
is  the  scattered  field.  Thus  g  is  known,  either  by  direct  measurement  of  Ey  on  the  interface, 
or  by  measurement  of  E*^^^  in  the  air  region  z  <  0.  In  Section  3  we  will  show  that  (3)  is  a 
strictly  hyperbolic  system  of  p.d.e. 
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The  incident  pulse  description  along  with  the  electromagnetic  boundary  conditions  on 
the  dielectric  interface  at  z  =  0  impose  symmetry  on  all  the  fields  and  a  change  of  dependent 
variables  (indicated  by  the  arrow)  is  possible:  Ey{x,  z,  t)  ->  Ey{z,  t-xjv),  and  i:,  t)  -> 

^  one-dimensional  system  for  the  tilded  fields  can  be  obtained  by  defining 
the  time-like  variable  ^  =  i  -  s/u,  and  changing  coordinates  (2,  z,  t)  -)■  (z,  in  (3).  Omiting 
details,  we  find  that  H,{z,0  =  Through  differentiation  with  respect  to  ^  and  z 

we  can  further  eliminate  to  obtain  the  following  one-dimensional  system  of  equations, 


d^Ey  _  , 

d^2 


dz^  €00  COS^  <f>inc  di 


M  pfl  pn 

E  " 


(4) 


dP^y 

di 


4.  ^  P"  —  E  •  n  =  1  M 


where  =  cf  yfiZ  cos  <f>inc  is  the  wavefront  speed  in  the  dielectric,  and.  c  -  1/ ^/e^.  In 
the  next  section  we  derive  from  (4)  a  single  equation  for  the  electric  field  in  the  dispersive 
medium  by  eUminating  all  the  Once  Ey  is  determined  the  remaining  magnetic  field  is 

Jo  Si  “s  • 

It  is  instructive  to  consider  here  various  limiting  values  of  the  r"  by  appealing  to  the 
solution  of  the  second  equation  in  (4),  i.e.. 


(0  =  ^-^Ey{^')<^'>  n  =  l, M.  (5) 

^ ^ 

Letting  some,  or  aU,  of  the  r"  be  large  corresponds  to  making  the  approximation  e"  /r"  ~ 
(1  _  0{^))It'^  in  (5).  By  letting  one  of  the  relaxation  times  be  large,  e.g.,  take  > 
r",  n  =  2, M,  so  that  (5)  becomes  Pj  ~  ^  Ey{(')dt\  one  can  model  the  presence  of  a 

zero-frequency  conductivity  in  the  medium  with  value  <r  =  The  right  hand  side 

of  (4)  will  have  the  term  -<r^  and  it  wiU  now  be  coupled  to  M  -  1  ordinary  differential 
equations  for  the  remaining  relaxations.  In  the  infinite  -relaxation  limit  (t  ^  l,7i  — 
so  that  P^  ~  Eyii'W)  (4)  becomes  a  telegraphers  equation  for  propagation - 

in  a  medium  of  constant  conductivity. 


d^Ey 


,  ^ 


+  O’c 


dEy 


0. 


(6) 


where  -  \ _  (E*Li  ^)-  In  this  case  the  “skin-depth”  extends  to  infinity.  Alter- 

€00  CO*  f^inc  fl  T 

natively,  one  can  derive  the  zero-relzccation  limit  by  integrating  (5)  by  parts  and  exploiting 
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the  fact  that  1/r”  -)■  oo;n  =  to  evaluate  the  Laplace  integral 

which  anises.  Then,  all  mechanisms  appear  in  equilibrium,  i.e.,  PJ*  ~  Ae"Py,  and  (4)  be¬ 
comes  a  lossless  wave  equation  describing  propagation  in  a  medium  of  nondispersive  relative 
permittivity  =  €«  cos’  (f>{ne  +  EiHi 

5’Py  _£L_^  =  0  17) 

de  i  + 

‘  too  eo»*  (pine 

Note  that  now  the  “skin-depth”  collapses  to  zero. 

Equations  (6)  and  (7)  Respectively  set  the  lower-  and  upper-bound  of  the  peak  electric 
field  inside  the  medium.  With  increasing  2,  Ey  becomes  smoother  since  and 

(6)  becomes  the  diffusion  equation  (no  propagation).  The  field  will  be  zero  past  some  small 
depth.  This  will  be  the  lower  bound  of  the  peak  response  in  the  original  medium.  On 
the  other  hand,  (7)  describes  the  upper  bound  of  the  response  in  the  equilibrated  medium. 
This  bound  is  the  zero  frequency  transmission  coefficient  for  the  interface,  T[uj  =  0)  = 
2/(1  +  v^).  We  have  just  shown  that  0  <  ^  T(a;  =  0),  i.e.,  that  the  maximum 

response  in  the  dispersive  medium  is  bounded  above  and  below  by  the  large-depth  response  in 
media  corresponding  respectively  to  the  “infinite”-  and  zero-relaxation  limits  of  the  original 
medium. 

3.  The  Wave  Hierarchy 

In  this  Section  we'  set  E  =  and  P”  =  P”  for  simplicity.  We  define  the  following 


operators: 


=  dt  +  — 


Dn  =  n  (®) 

m=l  * 
m^n 

M 

D  =  ■ 

m=l 

To  begin  the  derivation  of  (2)  apply  the  operator  D  to  both  sides  of  the  first  equation  in  (4) 
and  the  operator  to  the  second  equation  there  twice  differentiated  with  Our  working 


system  now  is 


5’P  2  -1 _ 
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n  =  1, M. 


Noting  from  (8)  that  D  =  I>„D"  we  sum  the  second  equation  in  (9)  over  n  thus  obtaining 
an  expression  for  the  right-hand  side  of  the  first  equation  there  in  terms  of  the  electric  field, 


M  ^2  pn  Af  1 

n=l  n=l 


Equation  (9)  now  becomes 


^  600  COS^  (pine  n=l  ^ 

Using  in  (11)  the  identity  ^Dn  =  (£>"  -  di)Dn  =  D-  lim.._oo{i:>"}^n  and  the  expansion 
^  =  E«naf-'‘ 


ao  =  1 


M  1  M 


t  M  -I  M  1  wi-lMI 

£>  =  —  — •••  E  ~  E  ”  =  i>  — >  M —1 

«•  ^  r’"  •  ^ .  t’""*  1  i -1 


I . .  In— . . . *l 


“Af  =  li 

we  obtain,  after  rearranging  the  various  terms,  the  sought  after  single  p.d.e.  for  E  in  the 
quarter-plcine: 

Af 

E  0.ef-'{d((  -  cld..)E  =  0;  z>0,(>0 


fio  —  1,  Co  —  Co 


1  M 

/3„  =  a„  H - r-T —  E  ~  ^  ~  •••>  ~  1 

Coo  COs2  0.^  ^  r*-»oo 


=  n  =  l,...,M-l 

VPn 

E^i  Ae‘  , 

“  '"V  /3m 
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The  signaling  problem  is  complete  once  we  specify  the  initial  and  boundary  conditions 
£;(z,0)  =  5{£?(z,0)  =  ...  =  d^'^^E{z,0)  =  0  and  =  g{()  respectively,  and  the 

radiation  condition  E{z  oo,^)  ->  0.  Equation  (13)  is  a  strictly  hyperbolic  p.d.e.  since 
the  n  =  0  term,  which  is  the  principal  pmt  of  the  operator  acting  on  E,  has  a  complete 
set  of  A/  +  2  distinct  eigenvectors,  one  for  each  of  the  two  distinct  eigenvalues  ico,  and  M 
eigenvectors  for  the  zero  eigenvzdue  of  multiplicity  M  due  to  the  operating  on  the  n  =  0 
term.  Note  that  (13)  also  holds  in  two-  and  three-dimensions  for  every  component  of  the 
electric  field  if  is  replaced  by  V^,  the  appropriate  higher  dimensional  Laplacian.  This 
is  because  the  incident  electric  field  induces  the  polarization  so  that  V  •  P  =  0  everywhere 
inside  the  medium.  Equation  (13)  is  vzJid  for  material  parameters  that  are  continuously  or 
discretely  layered  in  the  z-direction  (depth).  The  (3n  and  Cn  form  the  ordered  sequences 


/?0  =  1  </?!<...  <  /3n  < 

Co  =  V*’^**(cu  =  oo)  >  Cl  >  ...  >  Cn  >  ...  >  CM-\  >  Cm  =  =  0) 


(14) 


so  each  wave  type  is  a  member  of  a  wave  hierarchy. 


4.  Analysis  of  the  Hierarchy 

To  introduce  the  incident  pulse  duration  in  our  analysis  we  sccde  the  independent  Vciri- 
ables,  and  z‘  =  change  to  the  new  variables,  and  then  drop  the  primes  for 

convenience.  Our  working  equation,  factored  to  exhibit  left-  and  right-going  waves,  then 
becomes 

M 

'£,0.T;al‘-'{d(  -  i„a.m  +  5..a.)E  =  o,  (i5) 

n=0 

where  Co  =  1  and  c,»  =  c„/co.  The  “infinite”-relaxation  limit  (Equation  (6))  is  equivalent 
to  Tp  <  ^minimum  ^j^gnce  <  1>  while  the  zero-releixation  limit  (Equation  (7))  is 

equivalent  to  Tp  >  Tmaximum  >  1).  The  peak  response  will  be  bounded  as  described 

in  Section  2  for  any  pulse  duration  and  the  response  in  the  medium  will  depend  on  the 
frequency  content  of  the  boundary  function  g{().  The  highest  frequencies  constitute  the 
high-order  waves  described  mainly  by  n  =  0  term  in  (15)  and  will  propagate  with  the 
wavefront  speed  Co(=  1)  which  is  the  highest  speed.  In  the  “skin-depth”  the  high-order 
waves  approximately  satisfy  d^E  ~  — cod^E  since  the  relzixations  have  not  had  time  yet 
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to  respond  to  the  incoming  wave.  Using  this  in  each  term  in  (15)  except  in  d(  +  Cod,  we 
determine  the  p.d.e.  that  governs  this  order’s  propagation,  i.e., 


1  ^ 

(S(  +  lo8.)B  + 

^Po  n=l 


=  0, 


where  the  operator  is  interpreted  formally  as  n  -  1  integrations  of  E  with  respect  to  z. 

The  second  term  in  (16)  is  a  smoothing  operator  with  the  n  =  1  term  producing  exponential 
decay  of  E  with  depth.  Discontinuities  of  ^(0  will  be  important  only  near  the  interface 
across  the  characteristic  ray  z  -  co|  =  0  and  will  be  exponentially  attenuated  for  models 
with  M  =  1.  Models  fitted  to  data  with  M  >  1  will  smooth  any  discontinuity  in  the  incident 


pulse. 

Bsuids  of  lower  frequencies  2iround  the  relaxation  frequencies  are  each  governed  by  the 
intermediate  terms  1  <  n  <  M-1  in  (15).  Because  experimental  data  indicates  that  Pm  >  1 
{Pm  ~  0(10^^)  for  an  M  =  1  model  of  water)  these  intermediate  orders  are  unimportant  in 
comparison  to  the  low-order  waves  described  by  the  n  =  M  term.  From  the  dependence  of 
Pm  on  the  angle  of  incidence  we  also  deduce  that  the  ampUtude  of  the  low-order  waves  (and 
hence  of  the  main  response)  will  be  maximal  for  pulses  that  are  normally  incident.  Due  to 

the  largeness  of  Pm  the  low-order  waves  approximately  satisfy  d^E - and  using 

tHs  in  (15),  except  in  the  term  +  cmS,,  the  p.d.e.  governing  the  propagation  of  this  wave 


order 


(3«  +  c«8.)B+  1)  (j„)n  0. 


Equation  (17)  describes  diffusive-dispersive  propagation  with  speed  cm  (~  1/9  for  the  M  1 
water  model).  Dispersion  comes  from  those  terms  of  the  second  part  of  (17)  which  contain 
odd  number  of  derivatives  of  E  with  respect  to  z,  while  diffusion  comes  from  terms  with 
even  number  of  z-derivatives  of  E.  The  response  is  smooth  across  the  sub-characteristic  rays 


z  -  c„^  =  constant]  n  =  1, ....  M.  The  largeness  of  Pm  indicate?  that  low  frequencies  are 
important  in  such  problems.  These  frequencies  need  not  be  present  in  the  problem  through 
a  zero-frequency  component  in  the  incident  pulse,  rather  they  appear  through  the  spectrum 
of  the  envelope  of  incident  pulses  which  always  contains  low  frequencies.  One  could  have 
derived  (17)  by  using  the  slow  variables  z  =  (.''{z  —  c^/t)  and  t  =  (0  <  e  1)  in  the 

same  way  long-wave  approximations  are  derived.  In  such  a  derivation  the  particular  value 


of  7  will  depend  on  whether  M  -f  2  is  even  or  odd. 

In  this  paragraph  we  consider  (16)  and  (17)  for  M  =  1  (our  water  model).  For  this 
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medium  we  obtain  explicit  solutions  describing  the  behavior  of  the  response  and  we  are 
interested  in  deducing  qualitative  results  from  them.  These  results  will  be  vahdated  with 
a  numerical  simulation  in  the  next  section  for  4>xnc  =  0  (^  =  t).  The  equation  for  the 
smadl-depth  response  is 

(3;  +  d.)E  +  -  c\)B  =  0.  (18) 

while  for  large-depth  it  is 

(at  +  c-,a,)E  -  =  0.  (19) 

The  constants  zire  ci  =  ~  r(^  coA ~  ^oo-  Values  for  c,,  Coo, 

and  T  for  water  are  given  in  the  introduction.  Note  that  (19)  is  parabolic  while  the  (18)  and 
the  full  problem  are  hyperbolic.  The  solution  of  (18),  subject  to  the  boundary  condition 
-5(0,0  =5(0.  is 

E[z,i)  =  g[i  -  •  (20) 


It  shows  that  the  contribution  of  the  exponentially  decaying  pzirt  of  the  response  will  be 
negligible  after  a  depth  O(cot)  m  (in  dimensionzd  variables)  which  we  name  the  “skin-depth.” 
The  solution  of  (19)  is 


E{^z 


AT, 


I 


h{K) 


27r(l-5i)  Jo  U'-/c)> 


(21) 


where  is  related  to  the  inverse  Fourier  transform  of  the  spectrum  of  g[^)  after  the 

highest  frequencies  have  decayed  in  the  ‘‘skin-depth.”  The  origin  of  ^  is  the  first  time  after 
which  the  response  evolves  according  to  (21)  and  is  approximately  0(r)  sec.  The  response  is 
dififusing  around  the  sub-characteristic  z  =  ci^  ,  and  its  peak  will  be  on  that  ray.  Thus  the 
frequency  content  of  the  response  is  constricted  through  diffusion.  The  kernel  in  (21)  acts  like 
a  delta-function  because  ^  1.  This  leads  to  the  prediction  that  on  the  sub-characteristic 
ray  the  peak  response  for  an  arbitrary  incident  pulse  will  be  decaying  as  as  z  a  >  1/2,  or 
diS  ^  (this  can  be  obtained  by  remembering  that  also  ^  =  z/ci  on  the  sub-characteristic 
ray),  i.e., 

Ei 

For  pulses  with  /i(0)  =  0  (21)  is  evaluated  after  /i(/c)  is  expanded  in  a  Taylor  series  around 
zero. 


/ 


PiT,> 


.MO), 


27r(l-af) 


(22) 
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5.  Numerical  Experiments 

Now  we  validate  the  results  of  the  previous  section  for  the  M  =  l  Debye  model  of  water, 
namely  the  predictions  that  the  response  is  exponentially  decaying  in  a  thin  layer  (the  “skin- 
depth")  where  it  travels  with  the  infinite-frequency  phase  velocity  co  =  c,  and  that  past 
this  layer  it  travels  with  the  zero-frequency  phase  velocity  and  is  a  diffusion  wave  whose 
peak  decays  algebraically.  The  propagation  of  a  normally  incident  =  0°)  square  pulse 
of  Tp  =  50  psec  is  simulated  by  solving  the  system  (3)  with  M  =  1.  The  discretization 
parameters  are  chosen  witfr  the  aid  of  our  previous  work  [9]  and  the  results  shown  contain 
no  numerical  artifacts.  The  numerical  source  is  in  free-space,  5  cells  in  front  of  the  interface, 
and  there  are  10  cells  between  the  left-hand  side  radiation  condition  (which  is  in  free-space) 
and  the  dielectric.  We  record  time-traces  of  the  response  every  ten  spatial  cells  starting  at  a 
depth  of  10  cells  in  the  medium,  and  subsequently  search  in  the  data  for  the  peak  response 
which  we  save  along  with  the  location  of  its  occurence.  The  computational 

domain  is  taken  to  be  large  so  that  artificial  reflections  from  the  right-hand  end  of  the  grid 
do  not  reach  the  recording  locations  in  the  time  window  of  interest. 

Fig.  1  shows  the  (f ,  z)  location  of  the  peak  of  the  response  in  the  half  space  for  the 
first  80  psec  of  evolution.  The  z-axis  offset  at  f  =  0  corresponds  to  twenty  ceUs  (left-hand 
radiation  condition  to  first  recording  depth).  We  see  that  the  peak  response  indeed  travels 
with  the  free-space  speed  of  light  as  predicted.  The  prediction  is  graphed  as  a  dotted  fine. 
The  gap  in  the  curve  (between  the  parallel  Unes)  is  due  to  the  movement  of  the  peak  from 
the  leading  edge  of  the  pulse  to  the  trailing  edge.  The  data  in  that  interval  was  deleted 
since  the  apparent  slow-down  is  not  relevant  to  the  speed  of  propagation  of  the  peak.  Past 
50  psec  the  pulse  peak  starts  slowing  down.  This  happens  up  to  about  75  psec  whence  the 
speed  of  the  peak  has  achieved  the  zero-frequency  phase  velocity.  This  is  shown  clearly  in 
Fig.  2  where  the  evolution  of  the  peak’s  location  past  80  psec  is  shown.  The  slope  of  the 
stradght  line  is  c/\/80.35  and  this  is  very  close  to  the  observed  slope.  It  takes  a  “skin-depth” 
of  0{ct)  m  2.4  mm)  to  achieve  this  slow  speed. 

Fig.  3  shows  the  decay  ofthe  peak  of  the  response  for  0.05  <  z  <  0.25  mm  during 
the  early  evolution.  The  Figure  shows  that  the  predicted  exponential  decay,  i.e.,  the  dotted 
line  drawn  according  to  (20),  is  observed  in  the  simulation  results  (circles).  Finally,  Fig.  4 
depicts  the  decay  of  the  pulse  peak  for  z  >  0(cr)  m,  i.e.,  past  the  “skin-depth”.  On  the 
graph  we  indicate  with  a  dashed  line  the  algebraic  decay  z  5  predicted  by  (22)  with  ~  c^t. 
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X  Simulation 
-  Eq,  (18) 
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Figure  1  The  (t,  z)  location  of  the  peak  of  the  response  during  the  first  80  psec  of  evolution. 
The  slope  of  the  dotted  lines  is  exactly  c  as  predicted  by  Eq.  (18).  The  separation  of 
the  two  parallel  lines  is  exactly  the  pulse  duration. 
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Figure  2  The  {t,z)  location  of  the  peak  of  the  response  for  t  >  80  psec.  The  slope  of  the 
dashed  line  is  the  prediction  C\  =  c]  V80.35. 
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Figure  3  The  small-depth  decay  of  the  peak  of  the  response  as'  a  function  of  z.  The  dotted 
line  indicates  the  predicted  exponential  decay  and  the  circles  indicate  the  observed 
decay. 
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Figure  4  The  large-depth  decay  of  the  peak  of  the  response  as  a  function  of  z.  The  dashed 
line  indicates  the  predicted  algebraic  decay  and  the  circles  indicate  the  observed  decay. 
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The  observed  decay  is  plotted  with  open  circles. 


6.  Ramifications  for  Numerical  Schemes 

In  summary,  we  have  found  that  the  pulse  response  in  a  relaxing  dielectric  is  mainly  a 
diffusion  wave  traveling  with  the  zero-frequency  phase  velocity,  cm,  of  harmonic  waves,  and 
that  the  fastest  speed  co,  i.e.,  the  infinite-frequency  phase  velocity,  is  important  only  in  a 
thin  layer  near  the  air/di  electric  interface  in  which  the  response  is  hyperbolic  and  decays 
exponentially.  As  an  example,  for  the  M  =  1  model  of  water  permittivity  this  thin  layer  is 
0(10“^)  m,  and  ci  ~  co/9.  Thus,  one  is  faced  with  a  stiff  problem  in  the  time  direction  due 
to  the  exponentizJ  decay  in  the  “skin-depth,”  and  with  the  existence  of  dispeirate  wave  speeds 
in  well  defined  spatio-temporzd  regions.  In  addition,  the  problem  is  asymptotically  singidmly 
perturbed  since  for  large  z  and  t  it  changes  type  from  hyperbolic  (18)  to, parabolic  (19).  In 
this  section  we  indicate  how  these  findings  can  be  used  to  correctly  set  the  discretization 
in  existing  numerical  methods  for  pulse  propagation  in  relzixing  media.  Also,  we  present 
am  airgument  in  favor  of  using  high-order  finite  difference  schemes  which  aire  second-order 
accurate  in  time  and  fourth-order  accurate  in  space. 

We  have  shown  in  [9]  that  the  timestep,  At,  for  Leapfrog  based  Debye  schemes  is  required 
to  resolve  the  shortest  relaxation  timescaJe  r'”*”  for  reasonable  accuracy  over  long-time  sim¬ 
ulations.  This  is  due  to  the  time  stiffness  in  the  problem  for  realistic  media  and  to  the  fact 
that  simple  A-stable  approaches  (the  trapezoidal  method)  were  used  in  those  methods  to 
discretize  the  ordinary  differentiaJ  equations  for  the  polzurization. 

Next  we  show  how  to  chose  the  spatial  resolution  in  light  of  the  existence  of  multiple 
speeds  in  the  dielectric.  In  explicit  finite  difference  schemes  for  equations  with  multiple  wave 
speeds  the  stability  requirement  is  based  on  the  highest  speed  in  the  problem,  Co,  while  the 
spatial  resolution  is  based  on  the  slow  speed,  cm-  Due  to  the  large  magnitude  of  Pm  the 
slow  speed  will  be  dominant  in  realistic  models  of  matericJs.  In  the  system  (3)  the  other 
speeds' c„,  1  <  n  <  M,  me  not  evident.  It  must  be  noted  here  that  the  slow  speed  is  not 
characteristic  and  that  it  develops  during  the  evolution.  This  is  contrary  to  what  happens  in 
elastic  wave  propagation  where  the  two  speeds  (compressional  and  shear)  are  characteristic 
speeds.  If  one  considers  a  fixed  frequency  /  then  we  see  that  the  wavelength  associated  with 
the  slow  speed  is  smaller  than  that  associated  with  the  fast  speed,  thus,  a  small  length  scale 
develops  in  time.  If  we  assume  that  the  M+ 1  wave  types  are  decoupled  then  the  timestep  for 
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each  type  is  determined  by  Af/Az„  =  i//c„,  n  =  0, ....  M ,  where  :/  is  the  Courant  number. 
Thus,  for  a  fixed  timestep  and  Courant  number,  the  spatial  steps  for  the  most  important 
waves  in  the  problem  are  related  as 

AzAf  =  — Azq.  (23) 

Co 

We  deduce  that  the  disparate  speeds  require  a  reduction  of  the  Courant  number  since  one 
would  have  to  use  Azm  from  (23)  in  the  stability  restriction  CoAi/Az,vf  =  i/  to  obtain  the 
timestep,  i.e.,  coAt/Azo  =  ucmIco  <  i',  where  u  is  the  value  we  would  have  used  if  we  only 
knew  of  the  chciracteristic  speed.  Typically,  0.1  <  caz/co  <  0.5. 

Finite  difference  schemes  which  are  second-order  accurate  in  both  time  and  space,  the 
so  called  (2  -  2)  schemes,  require  that  the  Courant  number  used  be  the  maximum  allowed 
for  stability  (i/  =  1  in  one  dimension)  in  order  for  them  to  introduce  the  least  phase  error. 
Here  we  have  shown  that  once  the  spatial  resolution  is  set  according  to  the  slow  speed  one 
cannot  get  away  from  having  to  reduce  the  Courant  number.  However,  (2  —  4)  schemes  that 
are  second-order  accurate  in  time  and  fourth-order  accurate  in  space  are  stable  for  u  <  6/7, 
operate  well  for  0.1  <  i/  <  0.5,  and  they  are  overall  fourth-order  accurate  for  At  ~  0(1)  Az^ 
This  last  relationship  between  the  time  and  the  spatial  steps  in  the  high-order  scheme  is  like 
the  diffusion  scaling  t  ~  0(l)z*  which  appears  in  our  problem  asymptotically.  In  any  case, 
due  to  the  time  stiffness  one  desires  to  use  a  small  time  step  by  reducing  the  Courant  number 
or  by  increasing  the  spatial  cell  size  and  this  can  be  done  with  (2  4)  schemes.  If- (2  —  2) 

schemes  are  used  the  time  stiffness  will  also  affect  the  Az  by  forcing  it  to  be  extremely  small. 
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PHASE  ERROR  CONTROL  FOR  FD-TD  METHODS  OF 
SECOND  AND  FOURTH  ORDER  ACCURACY 


« 


Abstract 

For  FD-TD  methods  we  determine  the  spatial  resolution  of  the  discretized  domain 
in  terms  of  the  total  computation  time  and  the  desired  phase  error.  It  is  shown  that 
the  spatial  step  should  vary  as  Ax  g  in  order  to  maintain  a  prescribed  phase 

error  level  throughout  the  computation  time  tci  where  s  (=2  or  4)  is  the  spatial 
order  of  accuracy  of  the  scheme  and  ^  is  a  geometric  factor.  Significantly,  we  show  that 
the  thumb  rule  of  using  10-20  points  per  wavelength  to  determine  the  spatial  cell  size  for 
the  standard  scheme  is  not  optimal.  Our  results  are  verified  by  numerical  simulations 
in  two  dimensions  with  the  Yee  scheme  and  a  new  4th-order  accurate  FD-TD  scheme. 
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1.  Introduction 


In  finite  element  methods  for  the  frequency-domain  Maxwell’s  equations  the  electrical 
size  of  the  computationzJ  domain  is  related  to  the  quantity  kA  =  2irlNppu,,  where  k  is  the 
wavenumber,  A  is  a  typiczd  spatial  cell  size,  and  Nppu,  is  the  points  per  wavelength.  In  these 
methods,  in  order  to  maintmn  a  constant  error  level,  Npj„u  should  increase  as  the  domciin’s 
electrical  size  increases  [l]-[3],  thus  higher  order  elements  should  be  used  on  electrically  large 
problems  since  they  require  a  smaller  Npp^  for  the  same  error  in  compeirison  to  the  standard 
elements.  Herein,  we  consider  time-domain  FD-TD  methods  and  derive  the  Nyp^  required 
so  that  only  a  prescribed  amount  of  phzise  error  accumulates  in  a  given  computation  time 
interval.  We  will  consider  the  stand2Lrd  Yee  scheme  [4],  herezdter  referred  to  as  the  (2,2) 
scheme,  and  a  4th-order  accurate  in  space  FD-TD  type  scheme  [8],  hereafter  named  the 
(2,4)  scheme  since  it  is  2nd-order  accurate  in  time. 

We  will  derive  and  numericcdly  validate  the  estimate 

Wppu,  ~  a(3,0)(— ]  •  (1.1) 

e^ 

for  determining  the  spaticd  discretization  in  two  dimensional  FD-TD  schemes  of  spatial  order 
of  accuracy  s  (=2  or  4  here).  In  (1.1)  0  is  the  angle  of  propagation  with  respect  to  a  grid  axis, 
P  is  the  number  of  periods  in  time  for  which  the  computation  will  proceed,  and  is  the 
maximum  phase  error  in  radians  that  will  be  allowed  to  accumulate  over  the  computation 
duration  for  the  highest  frequency  in  the  problem.  Prom  tc,  the  actual  computation  time,  and 
a;,,  the  highest  frequency  that  is  expected  to  be  present  in  the  calculation  it  is  P  =  tcU;,/27r. 
An  estimate  for  w,  can  be  obtained  by  considering  the  initial  pulse  frequency  content,  while 
and  are  defined  a  priori.  Our  derivation  will  start  from  the  semi-discrete  Maxwell’s 
equations  thus  our  estimate  will  be  formzilly  valid  in  the  case  v  where  u  =  cAt/A 

is  the  Courant  number  of  the  fully  discrete  scheme.  However,  we  will  show  in  Section  2 
that  the  dependence  of  the  phase  error  on  P  and  Nppyu  as  shown  in  (1.1)  also  holds  for  fuUy 
discrete  schemes.  The  implication  of  (1.1)  is  that  methods  of  higher  order  of  accuracy  in 
space  can  achieve  the  same  error  over  a  fixed  computation  time  interval  as  the  Yee  scheme 
but  with  a  larger  spatial  step  (less  memory)  and  less  amount  of  timesteps.  (1.1)  will  be 
verified  by  numericed  simulations  of  mode  propagation  in  a  metallic  waveguide  (a  truly  two 
dimensioned  test  case)  with  the  (2,2)  and  (2,4)  FD-TD  schemes. 

The  first  (2,4)  and  (4,4)  schemes  on  two  and  three  dimensional  staggered  meshes  appeared 
in  [5]  for  subsets  of  Maxwell’s  equations  in  the  time-domain.  Reference  [6]  presented  a 
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method  for  achieving  4th-order  time  accuracy  without  requiring  the  introduction  of  extra 
time  levels  of  storage.  In  [7]  a  two  dimensional  (2,4)  FD-TD  type  scheme  was  presented  for 
the  elastic  velocity-stress  vector  wave  equations.  References  [8]  and  [9]  also  present  4th-order 
schemes  for  Maxwell’s  equations. 

2.  Dispersion  Analysis 

Our  working  equations  are  the  two  dimensional  Maxwell  s  equations  for  Transverse  Mag¬ 
netic  (TM)  polarization  i»an  unbounded  domain, 

dHy  _  dE, 
dt  dx 

=  JJ±  (2.1) 

dt  dy  • 

dEr  dHy  dH. 

dt  dx  dy  ’ 

but  our  analysis  can  be  applied  to  Transverse  electric  (TE)  polarization,  and  to  the  general 
three  dimensional  equations.  (2.1)  results  from  the  dimensional  MaxweU’s  equations  by  scal¬ 
ing  the  electric  and  magnetic  fields  respectively  on  and  s/iiot  ^^d  subsequently  scaling 
space  and  time  respectively  on  the  speed  of  light  c  and  its  inverse  c  so  that  |A|,  the  magni¬ 
tude  of  the  wavevector  k,  2ind  lj  axe  synonymous  quantities.  An  exact  solution  of  (2.1)  corre¬ 
sponding  to  a  spatial  Fourier  component  is  the  vector  M  =  {Hy{x,y,i),  Hx{x,y,t),  Ei{x,y,t)) 
=  *^^®  superscript  T  denotes  transpose,  and  fc  =  tfc*  + 

jky  is  the  wavevector.  We  determine  the  hy[t),hx{t),ez{t)  by  substituting  in  (2.1)  and  solv¬ 
ing  the  resultant  Ordineiry  DifFerentiad  Equations  (ODE)  for  the  time  dependence  with  initial 
conditions  given  by  the  constant  vector  M(x,y,0).  The  resulting  plane  wave  solution  is  given 
by 

f  /ly(0)  ) 

M  =  I  hx{0)  I  (2.2) 

I  e,(0)  J 

In  (2.2)  we  have  used  k*  =  |/!|sin^,  ky  =  |k|cosfl,  with  9  being  the  angle  of  propagation 
with  respect  to  the  vertical  y-axis  of  the  grid.  The  analysis  now  proceeds  by  considering 
2nd-  amd  4th-order  accurate  finite  difference  approximations  for  dfdx  and  dfdy  in  (2.1). 

Representing  the  spatial  derivatives  by  half-ceU  centered  2nd-order  accurate  finite  differ¬ 
ences  (FD-TD),  denoting  the  resulting  approximate  solution  by  and  substituting  it  in 
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(2.1)  we  obtzdn  the  ODE  system 


dt 


e  .  =-(. 


l)er(0 


dt 

ffe“PP(t)  1  .  ..^  _ 1  ,  .-ieA 


(2.3) 


.iky  A 


dt 


=  ^(e*^  -  e-‘^)h7(t)  -  -  e-'^)h:^{t) 


with  A  being  the  cell  wid^i  (taken  to  be  the  same  in  the  x  and  y  directions),  where  space 
is  discretized  as  (x,y)  =  (mA,nA)  with  integer  (m,n).  Upon  elimination  of  the  magnetic 
components  from  (2.3)  we  get  a  second-order  ODE  for  the  time  evolution  of  the  electric  field. 
Solving  this  ODE  we  find  that  (using  A  =  jfr|^~ ) 


Et^{x,y,t)  =  e,(0)e 


i(*.*+*»v-* io,e)  ] 


(2.4) 


The  phase  error  is  defined  as  the  ordered  difference  between  the  exact  and  numerical 
phase,  =  ^^xact  -  ^numerical-  Subtracting  the  phase  of  (2.4)  from  that  of  the  electric  field 
component  in  (2.2)  we  find  the  phase  error  produced  by  the  spatial  finite  differencing  to  be 


,  TTsin^. 

P*  =  sm(-rj^ - ) 

iVppti; 

TT  cos 

Py  =  M-Jj - )• 

Thus,  grows  /inear/y  with  time,  t,  for  a  given  discretization  that  is  fixed  at  the  beginning  of 
the  simulation.  For  an  actual  computation  one  would  set  \k\  to  be  the  highest  wavenumber, 
in  the  problem  for  which  we  want  a  preset  accuracy  to  be  maintained  up  to  the  end 
of  the  time  stepping.  Prom  (2.5),  using  t  “  tc  and  the  Taylor  series  of  the  sin  function 
for  TvfNpput  small,  and  keeping  only  the  first  term  in  the  series  for  which  in  (2.5)  is  non¬ 
vanishing  we  obtain  the  approximate  dependence  of  Npp^  on  the  allowed  phase  error,  the 
total  computation  time,  and  the  angle  of  propagation,  as  follows  (using  =  27rP) 

Nppu,  ~  (i)3  7r5(sm‘‘0  +  co3*9)^ {—)'^ .  (2.6) 
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The  phase  error  of  the  fully  discrete  scheme  can  be  found  using  the  methods  in  [12].  In  the 
notation  used  here  it  is  given  by 

UTT 


e*  =  -2,rPH  -  I 

^  I/TT  ^ 


(2.7) 


where  i/  is  the  Courant  number,  and  p.,  p„  are  as  in  (2.5).  It  is  immediately  obvious  that 
(2.7)  reduces  to  (2.5)  for  small  u  regardless  of  the  size  of  irfNpj^.  (2.6)  is  subsequently 
recovered  for  small.  Additionally,  for  smaU  and  j/  arbitrary  it  is  easy  to  show 

that 


N, 


ppw 


l/TT  \  ^  V  Q  iVpjnii 

so,  by  neglecting  4th  order  and  higher  terms,  (2.8)  reduces  (2.7)  to  (2.6)  regardless  of  the 
Courant  number  i/.  Thus,  (2.6)  also  holds  for  the  fully  discrete  scheme  with  arbitrary  i/  as 

long  as  xINpptu  is  smaill. 

The  O(A^)  accurate  half-cell  centered  discretization  for  the  spatial  derivatives  in  (2.1)  is 
given  by 

^le=(m+i)A  =  ■^lg(/m+l-/m)-^(/m+2-/m-l)l  (2-9) 

where  ^  =  x  or  y,  and  for  electric  variables  m  is  an  integer  while  for  magnetic  variables 
m  =  m  with  m  an  integer.  When  discretized  in  time  (2.1)  will  result  in  the  (2,4) 
scheme  given  by  the  difference  equations  [2]  below: 

J.i)  =  +  J.j)  +  -r[C(B?(i  +  i.i)  -  E;(i,}))  - 

+ 1)  =  nTh'.j + 1)  -  + 1)  -  sr(i,;))  - 


+  2)  -  £?(<,  j  -  1))1 


(2.10) 
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where  7  =  At/A,  (  =  9/8,  and  rj  =  1/24.  Altering  77,  (  one  obtains  the  (4,4)  FD-TD  scheme 
[8]  which  czin  be  analyzed  with  the  methods  herein.  The  semi-discrete  case  is  recovered  by 
letting  At  ->  0  in  (2.10).  Proceeding  as  for  the  (2,2)  scheme,  we  obtain  the  phase  error  for 
the  (2,4)  scheme 

«  TT 


9  .  ,  TT  sin  5 .  1  .  ,  Stt  sin  0 

p*  =  - ) 


8 


Py 


9  •  ( 


ppw 

TT  cos  9 


24 


1 


) - sin( 

/  0/1  ' 


Nppyu 

3:r  cos  6 


(2.11) 


24  "'  '  N„ 


). 


^pptu 

where  |jfe|  =  |/:,|  (the  magnitude  of  the  maximum  wavenumber  in  the  problem).  Setting 
[fc.ltc  =  27rP,  and  expanding  the  sin  functions  for  smeJl  vedues  of  irJNpp^  (2.11)  becomes 


(2.12) 


The  phase  error  of  the  fully  discrete  scheme  (2.10)  can  also  be  found.  In  the  notation 
used  here  that  phase  error  is  again  (2.7)  but  with  p*  and  py  now  given  by  those  in  (2.11), 
Arguments  similar  to  those  after  (2.7)  apply  here  too  for  the  relation  of  the  exact  error  to 
(2.12). 

Relations  (2.6)  and  (2.12)  are  the  main  result  of  this  paper  since  they  give  an  estimate  for 
Npput  based  on  a  priori  selected  computation  parameters.  Figure  1  shows  the  Npp^  required 
to  maintain  a  phase  error  of  0.1  radians  (~  5.73“*)  over  P  time  periods  of  computation 
for  a  harmonic  wave  propagating  along  the  grid  axis  (0  =  0°)  and  eJong  the  grid  diagonal 
{9  =  45").  The  benefit  of  a  high  order  spatial  differencing  scheme  is  apparent  from  this 
Figure  pcirticularly  for  long  time  computations  (P  large)  with  such  a  severe  phase  error 
restriction.  Significantly,  Figure  1  indicates  that  for  the  standard  Yee  scheme  the  thumb 
rule  of  Npjno  =  10  —  20  will  be  good  for  computations  up  to  only  P  =  3  for  a  phase  error  of 
5.73".  Figure  1  also  appbes  to  the  fully  discrete  schemes  used  with  any  u  whenever  irlNj^ 
is  smzill.  Since  in  our  ansJysis  i/  is  the  same  (~  0)  for  both  the  (2,2)  and  (2,4)  methods 
another  deduction  from  the  Figure  is  that  the  higher  order  methods  wiU  allow  a  Icirger  time 
step. 
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3.  Numerical  Validation 


We  consider  a  mode  propagating  in  a  metallic  air-fiUed  waveguide.  The  system  (2.1)  is 
solved  numerically  on  the  domain  0<z<l,  0<y<l,  0  <  t  <  tc,  together  with  Dirichlet 
boundziry  conditions  on  the  waveguide  wadis  Ez{x,y  =  0,  t)  =  E^[x,y  =  l,t)  =  0,  Hy{x,y  = 
0,t)  =  Hy{x,y  =  l,t)  =  0  for  0  <  X  <  1,  and  periodic  boundairy  conditions  adong  the 
direction  of  propagation  Ei{x  =  0,y,t)  =  Ez{x  —  l,y,t),  Hy{x  =  0,y,t)  =  Hy{x  =  l,y,t) 
for  0  <  y  <  1.  An  exact  solution  representing  the  th  mode  is  easy  to  find  with  standaurd 
methods  [11].  In  the  following,  is  the  longitudinad  wavenumber,  and  I  = 

1,2, ...  is  the  waveguide  mode  index.  This  model  problem  provides  a  truly  two  dimensionad 
test  of  our  ainadysis.  In  a  well  resolved  computation  the  truncation  error  is  negligible,  the 
only  error  soxirce  being  phase  misrepresentation  due  to  the  differencing.  Figure  2  shows  for 
the  Yee  scheme  the  effect  of  phaise  error  (only  5.12")  on  the  relative  error  for  a  computation 
in  the  waveguide  designed  with  the  anadysis  herein.  Although  visually  (Figure  2a)  the  phase 
error  seems  to  be  responsible  for  just  a  small  shift  in  the  computed  field  we  note  (Figure 
2b)  its  lairge  effect  on  the  relative  error.  The  reliability  of  interior  field  calculations  depends 
on  the  relative  error  they  contadn.  The  maiin  problem  here  is  that  the  zeros  of  the  field  au:e 
not  computed  correctly  due  to  the  (small)  phase  error  adthough  the  absolute  field  vadues  axe 
very  accurate. 

The  Dirichlet  conditions  for  the  standard  FD-TD  posed  no  problem  since  it  is  enough  to 
apply  the  boundary  condition  .B,  =  0  the  waveguide  wadis  by  prescribing  the  E^  nodes  at 
y  =  0,  1  and  0  <  x  <  1,  for  the  boundairy  condition  to  be  automaticadly  satisfied  to  2nd-order 
accuracy  by  the  differencing  strategy.  However,  the  implementation  of  boundairy  conditions 
for  the  (2,4)  scheme  is  not  so  stradghtforwaird.  The  spatiad  stencil  of  the  4th-order  method  is 
twice  ais  long  as  that  of  the  standaird  FD-TD,  thus  special  treatment  is  required  for  electric 
and  magnetic  nodes  adjacent  to  the  metadlic  waveguide  wadis.  After  alot  of  experimentation 
we  used  the  symmetry  properties  of  the  electric  and  magnetic  fields  with  respect  to  the 
waveguide  walls  (the  electric  field  is  an  odd  function  of  y,  while  the  magnetic  field  is  an  even 
function  of  y  with  respect  to  the  walls)  to  implement  the  Dirichlet  condition  for  the  (2,4) 
method  with  4th-order  accuracy.  Matching  the  Yee  scheme  2nd-order  differencing  to  the 
(2,4)  interior  stencil  in  order  to  compute  the  nodes  adjacent  to  the  computational  boundary 
proved  to  be  unstable. 

The  model  problem  is  discretized  with  the  (2,2)  and  (2,4)  schemes  in  order  to  obtaiin  the 
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Figure  2  a)  Electric  field  along  the  waveguide  at  P  =  10,  y  =  0.5,  computed  by  the 
standard  Yee  scheme  with  5.12'’  phase  error  (solid  Une).  b)  Error  relative  to  the  exact 
solution  solely  due  to  the  “small”  phase  error  allowed. 


Figure  2  a)  Electric  field  along  the  waveguide  at  P  =  10,  y  =  0.5,  computed  by  the 
standard  Yee  scheme  with  5.12"  phase  error  (solid  line),  b)  Error  relative  to  the  exact 
solution  solely  due  to  the  “small”  phase  error  allowed. 


fields  at  the  final  timestep  corresponding  to  time  tg  =  where  P  is  the  number  of  domain 

traversals  (periods)  of  the  mode  whose  phase  speed  is  The  total  number  of  timesteps 
is  NMAX  =  tc/At  =  P  *  Njyjnu  *  We  set  k  =  v/Stt.  Both  the  standard  FD-TD  and 
the  (2,4)  schemes  are  initialized  by  prescribing  =  0)  and  H^y[x^y^t  =  “Af/2) 

from  the  exact  solution  thus  obtaining  the  first  waveguide  mode  traveling  to  the  right  along 
the  x-axis.  The  numerically  computed  fields  should  exactly  reproduce  the  initial  conditions 
dStcT  the  completion  of  NMAX  timesteps  if  there  is  no  phase  error.  However,  phase  error 
accumulates  during  the  actual  computation,  and  we  measure  it  by  assigning  360/iVpp,y  degrees 
in  each  spatial  cell  and  using  linear  interpolation  to  find  the  field  node  at  time  t  —  which 
now  has  shifted  due  to  the  phase  error.  For  both  schemes  u  =  0.01  (<^  1/ \/2).  The  ceU  size 
is  determined  from  the  9  =  0''  curve  in  Figure  1  for  P  =  10  and  =  5.73^"  (=  0.1  radians), 
and  corresponds  to  =  32  for  the  (2,2)  scheme  and  =  8  for  the  (2,4)  scheme.  Since 
At  =  J^A  the  4th-order  scheme  required  one  fourth  the  cimount  of  timesteps  required  by  the 
standzurd  scheme  per  computational  period  (715.5  timesteps  for  the  (2,4)  scheme  hence  we 
will  measure  at  even  multiples  of  the  period). 

Figure  3  shows  for  the  Yee  scheme  the  predicted  (dashed  line),  obtained  from  (2.6)  with 
Q  =  0®,  against  the  c^Jculated  (plus)  phase  errors  as  functions  of  computation  time  measured 
in  integer  multiples  of  the  period  corresponding  to  the  ^  =  1  mode.  We  see  the  predicted  error 
was  slightly  Izu'ger  than  the  one  encountered  in  the  computation.  This  can  be  explained  by 
considering  that  the  particular  waveguide  mode  can  be  thought  to  be  composed  of  two  plane 
waves  traveling  to  the  right  at  an  angle  =  90°  —  tan  ^(1/2)  with  the  y-axis.  Therefore 
the  Nppu,  used  in  this  experiment  was  generous.  The  measured  phase  error  for  the  4th-order 
scheme  (stars)  can  also  be  seen  to  follow  the  predicted  phase  error  (solid  line)  obtained  from 
(2.12)  with  9  =  0°.  Again,  the  measured  error  was  slightly  less  than  the  prediction  for  the 
same  reason  as  for  the  Yee  scheme. 

4.  Summary 

In  this  paper  we  considered  the  phase  error  due  to  the  spatial  discretization  in  the  two 
dimensional  FD-TD  method  for  Maxwell’s  equations.  We  derived  how  the  cell  size  should  be 
chosen  so  that  only  a  preset  phase  error  accumulates  during  a  given  amount  of  computation 
time.  Truly  two  dimensional  numerical  simulations  validated  the  analytical  results.  Although 
the  derivation  is  formally  valid  for  cases  where  the  numerical  computations  are  performed 
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Figure  3  Computed  versus  P  for  the  Yee  scheme  (stars)  and  the  (2,4)  scheme  (plus), 
where  t/  =  0.01  in  both.  The  theoretical  curves  are  for  the  Yee  scheme  from  Equation 
(2.6)  for  =  32  (dashed  line  given  by  [e^l  ~  0.5783P),  and  for  the  (2,4)  scheme 
from  Equation  (2.12)  for  =  8  (solid  line  given  by  |e^|  ~  0.642P),  both  for  9  =  0®. 
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with  V  smaU,  we  have  shown  the  scaling  of  the  error  to  hold  for  arbitrary  i/  provided  Tr/i^- 
is  smaU.  Employing  a  small  i/  in  the  (2,4)  FD-TD  is  necessary  if  any  benefit  is  to  be  gained 
since  the  truncation  error  is  0{Ae  +  A^)  so  it  should  be  At  ~  A^  for  an  overall  4th-order 
accuracy  to  be  obtciined.  Also,  using  a  smaiU  u  will  allow  for  a  manageable  number  of  spatial 
grid  cells  in  the  FD-TD  for  very  high  frequency  applications,  for  dispersive  media  (see  [13]), 
and  for  any  problems  that  require  an  extremely  small  timestep.  It  has  been  determined 
that  the  number  of  spatial  ceUs  can  be  further  reduced  by  using  a  4th-order  accurate  spatial 

scheme  thus  allowing  electrically  larger  problems  to  be  addressed  given  a  fixed  quantity  of 

# 

computational  resources  and  of  permitted  phase  error. 
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A  General  Description  of  Propagation  and  Scattering  for 
Electromagnetic  Pulses  in  Dispersive  Media 


Thomas  M.  Roberts  and  Peter  G.  Petropoulos 


We  develop  some  new  methods  for  describing  pulse  propagation  for  general  disper¬ 
sive  media,  using  a  Debye  model  for  water  as  an  example.  Short-pulse,  long-pulse, 
short-time,  and  long-time  approximations  are  presented.  We  explain,  a  factor-of- 
nine  effect  in  the  speed  of  waves  in  water,  which  seems  to  have  been  previously 
unnoticed.  We  also  study  the  following  problem:  Knowing  only  the  peak  amplitude 
and  power  density  of  an  incident  pulse,  what  can  be  said  about  the  peak  amplitude 
of  the  propagated  pulse?  We  provide  sharp  upper  bounds  for  the  propagated  am¬ 
plitude  and  reduce  the  computation  of  those  bounds  to  a  calculator  exercise.  These 
bounds  may  be  useful  in  controlling  the  electromagnetic  interference  or  damage 
produced  in  dispersive  media. 
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1.  INTRODUCTION 

The  computation  of  electromagnetic  pulses  in  dispersive  media  is  a  highly  devel¬ 
oped  field.  For  inst2Lnce,  a  single  paper, ^  published  in  1976,  contains  numerics  for 
the  propagation  of  TE-  and  TM-polarized  electromagnetic  pulses  that  are  incident 
obliquely  on  an  inhomogeneous,  anomalously  dispersive  medium.  Computational 
electromagnetics  has  developed  so  extensively  since  1976  that  it  now  appears  that, 
given  enough  computer  resources,  one  can  compute  the  propagation  of  just  about 
any  single  pulse  through  just  about  any  single  medium.  But  these  studies  of  single 
pulses,  even  of  millions  of  single  pulses,  have  not  demonstrated  that  every  microwave 
pulse  travels  through  water^  with  one-ninth  the  speed  of  light  in  vacuum.  That  fun¬ 
damental  factor-of-nine  effect,  which  appears  to  have  been  unnoticed  until  now,  is 
established  here  by  studying  an  electromagnetic  wave  equation  and  its  scattering 
operators,  which  are  the  natural  places  to  find  broadly  applicable  rules  that  govern 
propagation.  We  will  show  that  the  time  of  arrival  of  transmitted  pulses  in  zinoma- 
lously  dispersive  media  is  related  to  a  slow  speed  given  by  the  DC  phase  velocity  in 
each  medium.  This  paper  has  several  other  new  results,  which  relate  to  the  widths 
and  peak  amplitudes  of  pulses,  and  to  quantities  that  resemble  power  density.  These 
new  results,  concerning  broad  classes  of  pulses,  are  vahdated  here  using  standard 
numerical  methods;  a  new  numerical  method  for  estimating  errors  is  also  developed 
and  used,  and  some  results  of  laboratory  experiments  on  pulse  propagation  in  a 
muscle- equivalent  material  are  explained.  Our  results  will  be  shown  to  be  helpful  in 
proposing  optimum  sample  lengths  to  be  used  in  Time  Domain  Spectroscopy  studies 
for  the  accurate  determination  of  the  infinite-frequency  and  static  permittivities  of 
Debye-type  dispersive  media.^“^  Also,  as  shown  in  Section  4,  our  results  can  form 
the  basis  for  a  sensitivity  analysis  of  the  dependence  of  the  medium  response  on 
the  parameters  obtmned  from  different  fits  to  the  same  band-limited  experimental 
data. 

This  work  was  done,  in  part,  to  assist  in  the  development  of  health-and- 
safety  regulations  for  electromagnetic  pulses  in  human  environments.  Our  goal 
was  to  develop  methods  to  support  the  regulation  of  basic  quantities  such  as  the 
peak  amplitude  and  the  power  density  of  incident  pulses,  so  that  it  would  not 
be  necessary  to  regulate  every  detail  of  a  pulse’s  time  trace.  Toward  that  end 
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we  formulated  a  problem  In  making  inferences  from  incomplete  data.  We  asked; 
Knowing  only  the  peak  amplitude  jmd  power  density  of  an  incident  pulse,  what  can 
be  said  about  the  peak  amplitude  of  the  propagated  pulse?  We  also  asked  what 
the  incomplete  data  would  imply  about  the  values  inside  the  dispersive  medium 
of  the  time  derivative  of  the  magnetic  field  H  and  of  a  different  quantity  that  is 
related  to  power  density.  Of  those  three  quantities— peak  amplitude  and  power 
density  and  dtH— it  is  dtH  whose  size  is  most  closely  linked  to  the  time  scale  of 
short-risetime  pulses;  further,  dtH  is  particularly  important  because  it  would  be 
largely  responsible  for  eleJtromotive-force  currents  in  any  circuit-like  structure  that 
is  inside  a  dispersive  medium.  Our  results  in  this  matter  of  incomplete  data  are 
quite  concrete.  We  will  show,  for  instance,  that  the  peak  ampUtude  of  the  electric- 
field  part  of  a  propagated  microwave  pulse  is  always  less  than  0.150  V/m,  for  depths 
greater  than  2.00  mm  in  water,  whenever  the  incident  electric  pulse’s  peak  amplitude 
is  less  than  1.00  V/m  and  its  power  density  is  less  than  5.29  x  10“^^Watt/m^, 
regardless  of  the  other  detaUs  of  the  microwave  pulse’s  time  trace.  We  have  similar 
results  for  power  density  and  for  dtH.  The  development  of  such  general  rules  for 
pulse  propagation  may  put  the  computational  basis  for  pulse-safety  standards  on  as 
firm  a  basis  as  for  the  existing  standards.®  for  continuous  waves  and  periodic  wave 
trains. 

This  paper’s  dispersion  models  and  time  scales  are  motivated  by  laboratory 
experiments.  We  use  a  one-term  Debye*  model  that  tits  laboratory  data  for  water 
up  to  100  GHz,  and  our  numerical  tests  involve  pulses  with  or  without  DC-frequency 
content  whose  time  scales  are  characteristic  of  short-pulse  radar.  Our  methods 
apply  to  all  other  Debye-like  media,  and  can  be  generalized  for  the  two-term  and 
five-term  Debye  models  that  fit  laboratory  data  for  muscle  and  muscle-equivalent 
materials.*”^  Our  methods  can  also  be  generalized  for  non-Debye  media. 


IV  -  3 


2.  FORMULATION  AND  RESULTS 


A.  PDEs 

The  equations  governing  the  scattering  and  propagation  of  an  obliquely  incident 
pulse  on  a  homogeneous  dispersive  half-space  occupying  z  >  0  axe  the  time-domain 
Maxwell’s  equations  for  the  fields  HxyHt^Ey.  This  set  of  equations  is  coupled  through 
a  polarization  current  {~^)  to  a  differential  equation  that  describes  the  evolution  of 
an  orientational  poleurization  (Py)  mechanism  of  Debye  type®  (a  releixation  process): 

+  Py  =  AeEyy  wherfi  Ae  =  e,  —  Coo,  St  and  Coo  are  respectively  the  zero-  and 
infinite-frequency  relative  permittivities,  and  r  is  the  dielectric  relaxation  time.  This 
o.d.e.  together  with  the  constitutive  law,  Dy  =  CoiSooEy  +  Py),  result  in  the  model 
frequency-domain  relative  permittivity  e(a;)  =  good-  where  Co  is  the  permittivity 
of  vacuum.  This  model  is  fitted  to  frequency-dommn  experimental  data  for  a  r2inge  of 
frequencies  uf  in  order  to  fix  the  various  medium  parameters.  Typical  values  for  water 
in  the  microwave  frequency  range  are  e,  =  80.35,  goo  =  1-00,  r  =  8.13  psec.  The 
phase  velocity  of  each  frequency  component  in  such  a  medium  is  — , 

with  c  being  the  speed  of  light  in  vacuum.  In  the  subsequent  amalysis  v*’^**(0)  and 
t;*’^**(oo)  will  arise.  FinaJly,  operational  considerations  fix  the  pulse  shape,  /,  and 
its  duration,  Tp. 

The  electric  field  incident  on  the  heilf-space  from  the  air  side  (z  <  0)  is  a 
plane  pulse  z,  t)  =  f{t  —  x  sin  (f>inclc~  z  cos  4>in-:lc)  of  duration  Tp.  We  eissume 

the  pulse  has  been  in  contact  with  the  interface  since  t  =  —  oo.  On  the  interface, 
z  =  0,  the  toted  electric  field  is  Py(x,  0,t)  =  g{t  —  xfv)  where  v  =  cf  sin^;„e;  the 
total  field  is  known  by  direct  measurement  of  either  the  field  on  the  interface  or  of 
the  scattered  field  in  z  <  0.  Defining  the  time-like  variable  ^  =  t  —  x/v  vre  find 

that  Py(x,Z,t)  =  J?y(0,  Z,^)  — >  Ey^Zf^)f  —  /fx,x(0>  ■Z,  ^)  ^*1 

that  Ht{z,^)  =  ij—Ey(z,^).  Changing  coordinates  (x,z,  t)  —>  (z,^)  in  the  resulting 
one-dimensional  system,  zind  eliminating  Hy  through  differentiation  with  respect  to  ^ 
zind  Py  by  using  the  operator  5^  +  we  obtain  a  single  third-order  pmtied  differential 
equation  for  the  electric  field  E  —  Ey  (shown  in  factored  form), 


di[d(_  -  c,dx){d^  +  c,d;)E  +  -  c,d,){d^  -f  c^d^^E  =  0,  z  >  0,  (2.1) 

T 

where  c„  =  c/(,y^cos^;„<,),  P  =  l-bAe/(goo  cos^  (}>inc),  and  Ci  =  Co/V^.  For  =  0, 


T^/  A 


(  =  U  The  signaling  problem  for  (2.1)  is  completed  by  giving  the  boundary  condition 
f;(0,^)  =  g{$),  and  the  initial  data  0)  =  0)  =  0)  =  0.  Hx  follows 

once  E  is  known,  and  H,  =  The  following  results  are  derived  in 

Subsection  3A, 

Equation  (2.1)  describes  the  propagation  of  all  possible  waves  of  different  or¬ 
ders,  and  their  corresponding  speeds,  that  can  be  excited  by  an  arbitrary  pulse.  The 
coefficients  exhibit  an  explicit  dependence  on  the  angle  of  incidence,  and  on  the  pa¬ 
rameters  that  describe  thj  medium.  It  is  a  strictly  hyperbolic^  partial  differential 
equation  since  the  principal  part  of  the  operator  has  real  distinct  eigenvalues  (three 
eigenvalues,  dcQ,  and  0),  and  a  complete  set  of  eigenvectors.  Causality  follows  from 
this  last  sentence.  The  characteristic  contributed  by  the  zero  eigenvalue  can  be  visu¬ 
alized  by  considering  that  d^  +  c^dx  =  when  cj  =  0.  The  main  feature  of  (2.1)  is  the 
two  wave  equations  exhibiting  distinct  speeds,  Co  and  ci.  Pulse  propagation  is  gov¬ 
erned  by  these  two  speeds  in  mutually  exclusive  spatial  regions.  Disturbances  mainly 
described  by  the  principal  part  of  the  operator  in  (2.1)  will  be  called  high-order  waves, 
while  those  described  by  the  remaining  operator  will  be  called  lower-order  waves.*  The 
speed  Cl,  while  not  a  characteristic  speed  (it  is  sub-characteristic,  ci  <  Co),  is  impor¬ 
tant  in  the  zmalysis  and  has  several  ramifications  for  experiments.  Also,  Ci  v  (0) 
and  Co  —  v*^**(oo);  i.e.,  the  main  disturbances  will  propagate  with  the  distinct  speeds 
which  are  equal  to  limiting  values  of  the  phase  velocity.  It  is  worthwhile  to  empha¬ 
size  that  experinaental  data  indicates  Ci.*^  Cq,  e.g.,  ci  ~  0.1116co  for  water  in  the 
microwave  range;  the  problem  is  stiff  so  the  pulse  travels  in  the  half-space  with  either 
of  two  speeds  that  Me  disparate. 

The  high-order  term  describes  the  dominant  behavior  for  depths  z  <  0[cot) 
m,  and  the  effect  of  the  lower-order  term  on  the  the  high-order  waves  is  an  exponential 
decay  with  z.  The  penetrating  pulse  propagates  with  speed  c.,  in  this  shallow  depth 
(~  10"*  m  for  water)  which  we  name  the  “skin-depth”  for  pulses  since  it  is  remi¬ 
niscent  of  the  weU  known  frequency-domain  concept.’  Prom  experimentally  obtained 
data  typical  of  tissue  r  =  O(10-‘’)  sec,  and  =  0(10).  Thus  f  is  large,  and  we 
expect  the  bulk  of  the  penetrated  pulse  to  travel  with  the  speed  cj  since  (2.1)  is  then 
approximately  -  c’E„  =  0.  The  main  disturbance  wiU  be  a  lower-order  wave. 
The  effect  of  the  high-order  term  on  the  lower-order  waves,  which  travel  with  speed 


Cl,  is  diffusive  in  character  and  importzuit  for  z  >  0{cor)  m.  The  main  response 
diffuses  around  the  ray  z,  =  Ci^  on  which  the  peak  of  the  response  is  found.  The 
peak  amplitude  on  the  sub-characteristic  ray  decays  as  1/v/^,  or  as  l/^/^  (for  fixed 
depth).  A  consequence  of  this  is  that  the  peak  of  the  energy-like  quantity  will 
decay  as  1/z  (1/0- 

The  response  will  also  depend  crucially  on  the  pulse  duration.  This  parameter 
appears  through  scaling  ^  with  Tp,  and  z  with  cTp.  Now  ^  and  Co  and  ci  are 

normalized  by  c.  Pulses  with  appreciable  zimplitude  most  often  have  Tp  10  ®  — 10 
sec,  so  ^  is  still  large.  Pulses  that  are  long  with  respect  to  the  relaxation  time 
^Tp  — >  00,  or  equivalently  if  r  — >  0)  will  propagate  unattenuated  in  the  half-space 
with  amplitude  equal  to  the  DC  value  of  the  frequency-domain  trmsmission  coefficient 
regardless  of  the  pulse’s  DC-frequency  content.  The  field  just  zifter  the  interface  (no 
“skin-depth”  since  CqT  -4-  0)  satisfies  a  lossless  wave  equation  with  speed  Ci.  On  the 
other  h£md,  very  short  pulses  (Tp  — >  0,  or  equivalently  if  r  — >  oo)  will  not  penetrate 
far.  In  this  case  the  electric  field  in  the  hatlf-space  (since  now  CqT  oo)  sees  a  high 
constant  conductivity  medium  thus  it  satisfies  a  telegraphers  wave  equation  whose 
far-held  is  the  diffusion  equation. 


B.  Green  Functions 

This  subsection  describes  some  rules  of  wave  propagation  that  are  derived  from 
time-domaiin  Green  functions.  The  history  of  these  Green  functions  is  reviewed  in 
Subsection  3B.  We  will  first  state  some  results  involving  upper  bounds  on  prop¬ 
agated  peak  amplitudes  and  power-density-type  quantities.  These  upper  bounds 
are  easily  computed  and  they  are  independent  of  the  detailed  nature  of  the  incident 
fields.  The  bounds  are  developed  for  normzd  incidence  here  and  for  oblique  incidence 
in  Subsection  3B.  The  present  section  concludes  with  a  description  of  wave  speeds 
and  with  brief-pulse  and  long-pulse  approximations.  These  rules  are  all  illustrated 
numerically  using  the  water  paxameters  e,  =  80.35,  £00  =  1.00,  and  r  =  8.13  psec 
from  Subsection  2A;  and  the  results  are  easily  generalized  to  other  Debye  media 
and  to  non-Debye  media.  The  necessary  derivations  and  numerical  validations  are 
in  Subsection  3B  and  Appendix  A. 

For  normal  incidence,  let  the  y-polarized  incident  electric  field  be  f[t  -  z/cq) 
in  the  air-filled  half-space  z  <  0.  In  the  water-filled  half-space  z  >  0,  the  resulting 
electric  field  is 


( - f(t--)+  r‘^"dsf{s)GMz,t  -  s). 

\6.15xl0~®m/  \  Co  J  Jo 

(2.2) 


F(z,t)  =  Ey{z,t)  =  exp 


The  Green  function  GE(z,t)  is  graphed  in  Fig.  1  for  several  depths  z  >  0.  For  the 
boundary  z  =  0,  Ref.  9  derives  GB(0,t)  =  exp[-t/(2.00  x  10“^^s)]/i[t/(2.05  x 

10"^’s)],  where  h  is  the  modified  Bessel  function  of  the  first  kind;  the  oblique- 

incidence  generalization  is  (3.12).  The  magnetic  field  ^fx(z,  t)  =  —  J^^^^dsdzEy(z,  s)/ [j.q 


also  has  a  Green-function  representation  similar  to  (2.2). 

Many  of  our  new  results  are  based  on  (2.2)  and  Fig.  1.  Safety  standards 
may  be  affected  by  this  type  of  analysis;  so,  in  Appendix  A,  we  show  how  to 
estimate  the  percentage  error  in  computations,  with  the  foUowing  results  for  our 
computations:  (1)  The  pointwise  numerical  error  in  Ge(0.500  mm,  t)  is  no  more 
than  1.70%  of  the  peak  value  (with  respect  to  t)  of  |Ge(0.500  mm,  t)l;  (2)  The 
pointwise  numericcJ  error  in  Ge(4.00  mm,  t)  is  no  more  than  0.800%  of  the  peak 
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value  of  1G'e(4.00  mm,  01;  and  (3)  For  intermediate  depths,  the  relative  error  in 
decreases  monotonically  from  1.70%  at  z  =  0.500  mm  to  0.800%  at  z  — 
4.00  mm.  The  closed-form  expression  for  Gb(0,0  is  exact.^ 

We  will  use  the  following  three  norms: 

|h(z,-)Ii=  /  dtlh(z,i)l 

Jo 

( (2-3)  • 

j/i(z,-)Ioo  =  least  upper  bound  of  l/i(z,-)l, 

where,  for  each  depth  z,  the  lezist  upper  bound  llh(z,*)ioo  is  evaluated  with  respect 
to  t.  Then,  for  each  depth  z  from  0.500  mm  through  4.00  mm. 


|£(..0I»  <  1/(01- «xp  (6.15  X  10-?^)  + 

regardless  of  the  detsuled  nature  of  the  incident  electric  field  f{t).  The  right  side 
of  (2.4)  is  easily  computed,  given  the  functions  i^i(2)  =  11Ge(^j’)12  and  F2{^)  — 
IGeC^.OIoo,  which  are  graphed  in  Fig.  2.  Inequality  (2.4)  defines  upper  bounds  on 
the  peak  amplitudes  of  E(z,  -).  This  inequality,  and  all  of  our  other  Green-function 
results,  are  validated  numerically  in  Subsection  3B.  The  upper  bound  on  the  right 
side  of  (2.4)  is  almost  attained  in  one  of  the  numerical  validations.  In  that  sense, 
the  upper  bound  is  sharp. 

We  will  now  show  how  relation  (2.4)  could  be  used  in  a  safety  standard  for 
the  peak  amplitudes  of  internal  electric  fields.  Suppose,  for  this  hypothetical  exam¬ 
ple,  that  it  has  been  determined  that  peak  electric  fields  must  be  no  greater  than 
0.200  V/m  at  depths  greater  than  1.00  mm,  and  no  greater  than  0.150  V/m  at 
depths  greater  than  3.00  mm.  That  hypothetical  internal-field  standard  is  trans¬ 
lated,  using  (2.4)  and  Fig.  2,  into  a  more  easily  regulated  standard  on  incident 
electric  fields  f{t).  The  easily  regulated  standard  is 


(0.202)  |/(•)l« 
F2{z)lf{-)h 


(2.4) 


T\/  -  Q 


(X  1 0'*  sec 
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0  L- ^ 
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0.0015 
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0.0035 
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Pig.  2.  The  L,  nonn  F,  (z)  and  the  L.  norm  F,  (z)  of  the  Green  function  Gb  (z.  •)  for 

water. .  There  norm,  and  E,n.  (2.4)  reduce  the  upper-bound  computation, 
to  a  calculator  exercise. 
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(1)  Peak  value  of  |/(01  —  0*740  V/m  or 

(2)  j^Peak  value  of  1/(01  <  40,  000  V/m  and  d0/(0l^  ^  4.00  x  10 

#  oo 

(3)  fpcak  iraJue  of  !/(()!  <  60,  000  v/m  and  /  d(|/(t)|  <  2.80  X  10"*‘Vs/m  . 

(2.5) 


or 


By  reading  the  two  graphs  in  Fig.  2  and  using  a  calculator,  one  can  use  (2.4)  to 
show  that  any  incident  field  that  satisfies  item  (1)  or  item  (2)  or  item  (3)  of  (2.5) 
is  guaranteed  to  produce  internal  fields  that  comply  with  the  hypothetical  internal- 
field  standard,  regardless  of  all  other  details  of  f{t). 

Upper  bounds  also  exist  for  quantities  that  resemble  power  densities.^°  In 
particular,  for  any  incident  field  f{t)  and  for  each  depth  2  from  0.500  mm  through 
4.00  mm,  the  power- density- type  term^°  [11-S(2)  =  Jq  dt\E{z,t)\^  satisfies 


/~dt|E(2,t)p 

Jo 


< 


(0.202)  1/(012,1 1" 

F.(^)l/(-)l.  Jj  ■ 


(2.6) 


Subsection  3B  numerically  validates  the  inequality  in  (2.6),  showing  that  the  upper 
bound  on  the  right  side  of  (2.6)  is  almost  attained  in  at  least  one  case. 

We  will  now  show  how  the  upper  bounds  in  (2.6)  could  be  used  in  a  safety 
standard  for  an  power-density-type  quantity  related  to  internal  electric  fields.  Sup¬ 
pose,  for  this  hypothetical  example,  that  it  has  been  determined  that  the  power- 
density-type  quantity  /g°°dtlFr(2,  t)p  must  be  no  greater  than  1.50  x  10  V  s/m 
at  depths  greater  than  1.00  mm,  and  no  greater  than  1.00  x  10  V^s/m^  at  depths 
greater  than  3.00  mm.  Equation  (2.6)  and  Fig.  2  translate  this  hypothetical  stan¬ 
dard  for  internal  fields  into  a  more  easily  regulated  standard  on  incident  electric 
fields  f{t): 


-  II 


(2.7) 


I  VC?  or 


(1)  /  dt\f{t)\^  <  2.40  X  10-“V^s/i 

Jo 

(2)  |^jf*dt|/(0r  <  36.0  Yh/n?  and  jf  dil/(i)|  <  3.90  x  10"^^  Vs/m  . 


Subsection  3B  shows  that  some  pulses  almost  attain  the  upper  bounds  In  (2.6), 
which  W21S  used  to  obtsun  (2.7). 

f 

The  upper-bound  concepts  in  (2.4)  and  (2.6)  are  easily  extended  to  mag¬ 
netic  fields  and  their  time  derivatives,  and  to  oblique  incidence.  Subsection  3B  has 
numerical  results  for  all  of  those  extensions.  One  can  see  there  that  making  the 
angle  of  incidence  more  oblique  wiU  decrease  the  penetration  into  the  medium  of 
power-density-type  quantities  and  also  peak  electric  and  magnetic  fields.  A  related 
closed-form,  modified-Bessel-function  expression  for  the  oblique-incidence  reflection 
kernel  R^{t)  =  GB,«(0,t)  is  given  in  (3.12). 

We  will  now  state  some  Green-fimction  results  concerning  wave  speeds. 
These  results  are  derived  in  Subsection  3B.  For  simplicity,  the  results  atre  stated 
for  normal  incidence.  Our  first  conclusion  is  that  the  main  bulk  of  an  electromag¬ 
netic  pulse  travels  through  water  with  speed  cq  for  0.3  mm,  and  then  slows  until,  for 
all  depths  beyond  0.7  mm,  the  pulse  travels  with  the  constant  speed  co/9.0.  That 
behavior  contrasts  with  the  wavefront  speed,  which  is  mathematically  well  defined 
but  is  not  always  observable  in  a  laboratory.  The  wavefront  speed  is  precisely  cq 
for  all  depths.^ ^  Section  4  discusses  various  Debye  models  that  are  consistent,  to 
within  about  10%,  with  the  bamd-limited  water  data^  used  here.  The  laxge-depth 
speeds  (all  «  co/9.0)  for  those  Debye  models  vary  by  only  about  10%.  The  shallow- 
depth  speed  and  the  wavefront  speed  of  any  Debye  model,  however,  axe  both  equal 
to  Co  =  l/V/^o^oo-  The  shallow-depth  and  wavefront  speeds,  consequently,  change 
considerably  as  one  varies  the  Debye  parameter  Coo  from  1.00  through  10.0,  as 
described  in  Section  4.  Therefore,  a  measurement  of  the  wavefront  speed  or  the 
shzJlow-depth  speed  would  determine  the  Debye  pcirameter  Coo- 

We  conclude  with  some  Green-function  results  concerning  short-pulse  and 
long-pulse  approximations.  The  propagation  of  any  finite-valued  incident  electric 
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pulse  /(f)  is  given  simply  by  (2.2)  and  Fig.  1.  We  get  additional  insight  by  con¬ 
sidering  approximations  for  what  we  will  call  elemental  pulses  /*:  An  element  /e(f) 
is  zero  except  on  a  single  time  interval,  during  which  it  is  either  strictly  negative 
or  strictly  positive.  For  instance,  a  square  pulse  is  an  element,  but  a  one-cycle 
sinusoid  is  not  an  element.  Elements  are  important  because  any  incident  pulse  / 
is  a  sum  of  positive-valued  elements  and  negative-valued  elements.  If  the  dura¬ 
tion  of  an  element  /.(f)  is  much  briefer  than  30  psec,  then  the  propagated  pulse 
element  is  apprbximately  J/e(’)li^B(^>f)  (s®®  Fig.  1)  for  all  depths  greater  than 
0.7  mm.  If  the  duration  of  an  element  /^(f)  is  much  longer  than  50  psec,  then  the 
propagated  pulse  element  is  approximately  0.2/e  [f  9.0(z  1  mm)/c  -f- 17  psec] 

0.2/e  (f  -  9.0z/c  -  13  psec)  for  aU  depths  from  1  mm  through  the  depth  at  which 
the  duration  of  Ge(z»<)  becomes  comparable  to  the  duration  of  /e(f)..  These  ap¬ 
proximations  are  derived  in  Subsection  3B. 


TA/  -  n 


3.  DERIVATIONS 


A.  PDEs 

To  extract  from  (2.1)  the  equation  describing  the  emly-time  evolution  (in  the  “skin- 
layer”)  of  the  response  we  set  everywhere  in  (2.1)  ~  —Cod^  except  in  the  operator 

d^  +  Codz  since  it  expresses  the  propagation  of  the  high-order  waves.  Any  other  terms 
in  the  resulting  equation  will  describe  the  effect  of  the  lower-order  waves.  The  main 
disturbzmce  for  emly  times  is  modeled  by 

(Sf  +  =  Oi  z<c.T,  (3.1) 

r  2c^ 

subject  to  the  boundary  condition  E(0,^)  =  5(^).  The  solution  of  (3.1)  is 


=  g((  -  (^)]  ■  •  (3-2) 

We  see  that  the  response  decays  exponentiadly  in  a  thin  region  of  depth  z  ~  0{cot), 
where  the  speed  of  propagation  is  c,,.  Note  that  the  decay  constant  is  inversely 
proportional  to  cos^  (pine  thus  normal  incidence  will  result  in  the  greatest  amplitude 
in  the  medium.  To  describe  the  evolution  of  the  lower-order  waves,  which  travel  with 
speed  Cl,  we  set  in  (2.1)  d(  ~  —cidi  except  in  the  operator  +  cidz  which  expresses 
the  hyperbolic  nature  of  the  lower-order  waves.  The  mmn  disturbzmce  is  now  modeled 
by 

la,  +  cA)B  =  ^^^a]B-,  z>c,r.  (3.3) 


{ac  +  cA)E  =  -z-^a^.B-,  z>c,r.  (3.3) 

The  boundary  condition  is  approximately  B{zo,  ^*)  =  where  Zg  is  the  depth  after 
which  (3.1)  no  longer  applies,  is  the  time  with  origin  at  Zo/co  (the  time  it  takes 
for  the  pulse  to  reach  Zg  in  the  “skin-depth”),  and  )  represents  (3.2)  evaluated 
at  Zg.  Equation  (3.3)  is  an  advection-diffusion  equation,  and  describes  the  response 
aifter  a  depth  of  O(cor)  m.  The  pezik  of  the  response  is  on  the  sub-chziracteristic 
ray  =  ci^*.  The  solution  is  very  easily  obtained  from  the  solution  of  the  diffusion 
equation.  It  is^’ 


fi'  h(/c) 

dn - -  -^-r  exp 

“  (r  -  «)5 


-  ci(^*  -  «)P 
e-K 
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Various  techniques  can  be  used  to  estimate  the  integral  in  (3.4)  since  PJr  is  large. 
Here  we  are  interested  only  in  the  primary  behavior  of  i?  as  a  function  of  depth.  For 
^  ^  1  the  response  is 


E‘(z,()  =  MO)^  -  cj)  f  I  {  [2r((/-c5).  f 


On  2.  =  Ci4'  we  find  that  max{E‘}  ~  1/v/i,  or  ~  1/^.  This  is  verified  with  the 
numerical  experiments  in  Section  4. 
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B.  Green  Functions 

This  subsection  derives  our  Green-function  results  in  the  order  in  which  they  are  de¬ 
scribed  in  Subsection  2B.  Equation  (2.2),  for  insteince,  is  a  Green-function  represen¬ 
tation.  These  Green  functions  have  become  a  standard  technique  in  computationcd 
electromagnetics.  They  were  first  developed  for  non-dispersive  media, and  were 
then  used  to  compute  fields  in  dispersive  media.^^  That  dispersive-medium  work 
has  not  yet  been  published,  owing  to  the  death  of  R.  Krueger,  but  generalizations 
eire  available.^®  The  Green  function  programs  used  in  the  present  paper  were 
developed  by  the  authors  of  Ref.  16,  along  the  lines  of  Appendix  A  of  that  paper. 
The  Loo  and  Lj  norms  in  (2.3)  have  special  physical  significamce.  The  Lj 

norm  is  important  because  (ceo/2)  (l£?|l2)^  is  the  power  density,  whose  mks  units 
are  Watt/m^,  of  an  electric  pulse  in  free  space.^°  Consequently,  we  will  focus  on  the 
peak-value  (p  =  oo)  and  power-density  (p  =  2)  cases  of  the  inequcdities 

|E(z,.)|,<e-“|/(-)|,  +  lG(z,.)|,l/(-)l,.  (3.6) 

which  axe  obtained  by  applying  the  Young  theorem^*  to  (2.2),  and  for  which  a  = 
1.63  X  10^  m"^.  In  particular,  although  (3.6)  is  valid  whenever  1  <  p,  g,  r  <  oo 
satisfy  r~^  =  1  ■+•  p~^  —  q~^ ,  we  are  most  interested  in  the  cases 

IIGb(^,*)IiII/(*)1~> 

5E(z,.)Joo  <e-“'l/(.)ioo  +  min  lGB(^,-)bl/(-)l2,  (3.7) 

L  I|Gb(2.,-)Bco!/(*)1i  J 

and 


iGB(^.-)bj/(-)«i  Jj  ■ 


(3.8) 


In  numerical  computations  for  water,  BGb(2,')Ii  observed  to  decrease  slowly 
and  monotonically  from  0.2019  at  0.5  mm  to  0.2014  at  3.5  mm.  It  is  as  if,  for  those 
depths  in  water,  the  advection-diffuslon  equation  (3.3)  were  approximated  by  a  heat 
equation  and  GE{z,t),  which  is  positive  valued  for  depths  beyond  0.5  mm,  were 
anedogous  to  a  temperature  distribution  whose  total  conserved  heat  is  proportional 
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to  |C?b(2,*)Ii-  Equations  (3.6)-(3.8),  the  almost-constant  nature  of  ^nd 

numerical  Green-function  computations  produced  the  results  in  (2.4)-(2.7). 

We  now  consider  five  numerical  examples  that  vahdate  the  inequalities  in 
(2.4).  We  will  see  that  the  minimum  upper  bounds  in  (2.4)  are  almost  attained  for 
some  inddent  pulses.  For  these  five  examples,  we  choose  the  following  hypothetical 
restrictions  on  the  incident  electric  pulseS  /(t): 

ll/(*)loo  <  1-00  V/m  and 

4 

!!/(•)  II2  <  6.32  X  10"*  Vs^/^/m  and  (3.9) 

|/(•)ll  <4.00  x  10"“  Vs/m. 

Then  (2.4)  shows  that  any  incident  pulse  /  that  satisfies  (3.9)  will. produce  an 
internal  field  whose  peak  amplitude  satisfies 


IjE?(z,0Ioo  <  expf - — — s— )  V/m  +min  (6.32  x  10~®  Vs^/^/m)  Fi(z), 

*  ^  -  L^V6.15  X  10-^mJ  J  (4.00  x  10"“  Vs/m)  ^^(z) 


0.202  V/m, 


(3.10) 


where  J^i(z)  and  1^2  graphed  in  Fig.  2.  Each  sum  of  the  exponential  in 

(3.10)  and  a  term  from  the  “min”  dause  of  (3.10)  yields  one  of  the  three  top-most, 
boldface  curves  in  Fig.  3.  Relation  (3.10)  guarantees  that  the  depth-dependent  peak 
amplitudes  of  J5(z,-)  are  less  than  the  minimum  of  the  three  boldface  upper-bound 
graphs.  That  prediction  was  tested  using  five  incident  pulses  /  that  comply  with 
(3.9).  Those  incident  pulses  are:  (1)  a  40-psec-duration  square  pulse  with  1-V/m 
amplitude;  (2)  the  absolute  value  /a  (t)  =  1/4(01  of  the  4-cycle,  80-GHz  sinusoid  in 
item  4  below;  (3)  the  absolute  vcdue  —  l/5(t)l  of  the  1-cycle,  80-GHz  sinusoid 
in  item  5  below;.  (4)  a  4-cycle  80-GHz  sinusoid  with  1-V/m  amplitude;  and  (5)  a 
1-cycle  80-GHz  sinusoid  with  1-V/m  amplitude.  The  norms  of  the  pulses,  which 


are  tabulated  below,  <ill  satisfy  (3.9). 
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Peak  Value  (Volt/meter) 


Fig.  3.  Five  numerical  validations  of  the  upper-bound  concept  for  the  depth-dependent 

peak  amplitudes  of  The  boldface  curves  are  upper  bounds  from  re¬ 

lation  (3.10). 


IV  -  18 


Table  1.  Five  incident  pulses  that  satisfy  the  conditions  in  (3.9).  The  first  incident 
pulse  is  1  V/m  for  0  <  t  <  40  psec,  and  it  is  0  for  all  other  times.  The  second  pulse 
is  the  absolute  value  of  (1  V/m)  sin[27rt/(8  x  10'°s)]  for  0  <  t  <  50  psec,  and  it  is  0 
for  all  other  times. 


Example  Duration 

Type 

11  •  llco 

11  •  112 

I -111 

(V/m) 

(Vs^/^/m) 

(Vs/m) 

1 

40-psec 

*  square  pulse 

1.00 

6.32x10"® 

4.00x10"^^ 

2 

4  cycles 

|80-GHz  sine| 

1.00 

5.00x10"® 

3.18x10"“ 

3 

1  cycle 

|80-GHz  sine| 

1.00 

2.50  X  10"® 

7.96x10"“ 

4 

4  cycles 

80- GHz  sine 

1.00 

5.00x10"® 

3.18x10"“ 

5 

1  cycle 

80-GHz  sine 

1.00 

2.50x10"® 

7.96x10"“ 

The  peak  amplitudes  of  the  five  internal  fields  E{z,  •),  corresponding  to  the  above- 
tabulated  incident  fields,  are  also  graphed  in  Fig.  3;  the  curves  for  Examples  4 
and  5  almost  overlap.  Those  peak  amplitudes  are  all  less  than  the  (boldface)  upper 
bounds  described  earlier.  The  five  examples,  therefore,  numerically  validate  the 
upper  bound  concept  in  (2.4).  Fig.  3  also  shows  that  the  upper  bounds  are  sharp 
in  the  sense  that  the  peak  amplitude  of  one  pulse  (Example  1)  almost  attains 
the  minimum  upper  bound.  That  example  involves  a  pulse  with  a  nonzero  DC- 
component^*»-2i  dt/(f).  It  makes  intuitive  sense  that  the  presence  of  a  DC 
component  in  a  pulse  would  tend  to  diminish  the  attenuation  of  the  pulse  in  any 
medium,  as  an  elementary  analysis^^  affirms  for  a  non-Debye  medium. 

Having  just  validated  the  upper-bound  concept  (2.4)  for  peaks,  we  now  val¬ 
idate  (2.6):  The  two  top-most,  boldface  curves  in  Fig.  4  correspond  to  the  upper 

bounds  (e-^i/lh  +  llGfEllil/iz)"  and  (e-^^ i/lU -b IGElbl/ii)'  in  (2.6),  subject  to  the 
hypothetical  restriction  (3.9).  The  other  five  curves  represent  the  power-density- 

type  quantities  [lFlr(2,  produced  by  to  the  five  incident  fields  in  Table  1.  We 


Time  Integral  (Volt^sec/m^) 


Fig.  4.  Five  numerical  validations  of  the  upper-bomid  concept  for  the  power-density- 

type  time  integral  /~dilF7(2,  t)p.  The  boldface  upper  bounds  are  from  rela¬ 
tion  (2.6). 
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see  that  the  field  produced  by  Example  1  of  Table  1  almost  attains  the  minimum 
upper  bound  in  Fig.  4.  This  completes  our  validation  of  the  upper-bound  concepts. 

As  explained  above  (2.1),  oblique  incidence  is  taken  into  account  using  a 
simple,  widely-known  transformation  of  vciriables.  Using  the  transformation,  we 
obt^n  numerical  results  for  45“-incident  electric  fields  f[t  -  (x  -f  z)/(-\/2c)].  In 
water  (z  ^  0)  the  y  component  of  the  electric  field  is 


E{x,2,t)-E(o,z,t 
E{Q,z,t)  = 


— z 


.35  X  lO-^m 


M_  \  ft-z/{V2e) 


(3.11) 


The  function  Gb,45»  {-z,  t)  is  graphed  in  Fig.  5  for  several  depths  z  >  0.  The  numeri¬ 
cal  errors  in  Fig.  5  were  quantitatively  estimated  using  the  method  of  Appendix  A. 
The  results  are:  (1)  For  each  depth  z,  from  0.160  mm  through  2.88  mm,  the  er¬ 
ror  in  the  computed  values  of  Ge,45*('Z>0  more  than  3.19%  of  the  peak 

value,  with  respect  to  time,  of  the  actual  values  |GB('Z,t)lj  and  (2)  The  relative 
errors  decrease,  but  not  necessarily  monotonically,  from  3.19%  at  0.160  mm  to 
2.07%  at  2.88  mm.  At  the  boundary  z  =  0  and  for  all  t  >  0,  GE,45'>(0,i)  = 
-t“^exp[-i/(1.01  X  10“^^s)]/i[t/(1.03  x  10"^^s)],  where  h  is  the  modified  Bessel 
function  of  the  first  kind.  More  generally,  the  oblique-incidence  transformation  and 
Ref.  9  imply,  for  all  t  >  0,  that 


R’{t)  =  Ge,.(0, ()  =  -^ «P [-  (<>  +  2eSco^)  ‘] 

That  exact  result  uses  the  modified  Bessel  function  of  the  first  kind  to  represent 
the  reflections  in  the  air-filled  half-space  (z  <  0)  that  are  caused  by  waves  that  are 
incident  obliquely  on  the  Debye  half-space  (z  >  0)  defined  by 


D{ 


£oo^0.®('Z,  t), 

£oo£oE{z,t)  -h  aeo/od3e"*’(‘“'^E(z,3), 


z  <  0 
z  >  0  ■ 


(3.13) 
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Fig.  5.  The  time  dependence  of  the  Green  function  (■*>0  several  depths  z. 


This  Green  function  is  used  for  pulses  that  are  incident  on  water  at  an  angle 
of  45*. 


In  particular,  the  reflected  field  is  exactly  ds/(s)i2®(t  -  s)  for  any  incident  field 
/.  Equation  (3.12)  was  used  numerically  to  validate  the  z  =  0  boundary  values  of 
the  oblique-incidence  Green-function  computations  that  yielded  Figs.  5,  6  and  8. 
Equation  (3.12)  was  also  used  to  validate,  at  z  =  0,  the  previous  percentage-error 
estimates  for  oblique  incidence.  At  z  =  0,  the  estimated  relative  error  was  1.91%; 
the  true  relative  error  was  1.87%. 

Applying  the  Young-theorem  result  (3.6)  to  the  oblique-incidence  represen¬ 
tation  (3.11)  yields  obvi(^us  oblique-incidence  generalizations  of  the  upper-bound 
results  (2.4)-(2.6);  for  instance,  |E(0,z,-)i2  <  |/(•)l2  exp[-z/(4.35  x  10-®m)]-|- 
rmii[iGB,45«(^,-)iil/(*)ll2»  iG'E,45»(2,-)l2l/(*)li]-  ThenormsE3(z)  =  11Ge,45*(^, Ob 
and  Fi{z)  =  iGB,45»(2,-)loo  are  graphed  in  Fig.  6,  and  iGB,45»  (-z,  Oil  was  observed 
to  decrease  slowly  and  monotonically  from  0.149  at  0.240  mm  to  0.144  at  3.00  mm. 
We  numerically  tested  these  oblique-incidence  inequalities  using  the  five  pulses  in 
Table  L  The  inequalities  were  vzJidated  in  each  case,  and  the  minimum  upper 
bound  was  almost  attained  in  the  case  of  a  DC-component  pulse  (Example  1). 

We  will  now  substantiate  the  results  in  the  last  two  paragraphs  of  Subsec¬ 
tion  2B.  The  results  rely  mainly  on  (2.2)  and  Fig.  1.  Note  that  the  first  term  on 
the  right  side  of  (2.2)  represents  a  wave  that  travels  with  speed  cq  and  decays  ex¬ 
ponentially  by  a  factor  of  132  in  each  0.300  mm  interval.  Therefore  the  convolution 
term  in  (2.2)  predominates  for  depths  greater  than  0  300  mm.  The  major  features, 
such  as  the  peaks,  of  the  convolution  kernel  Ge  ^.re  seen  in  Fig.  1  to  travel  more 
slowly  than  cq  for  depths  greater  than  0.500  mm.  The  peak  of  Ge(2^)0j 
stance,  is  shown  in  Fig.  7  to  travel  with  speed  cq  for  the  first  0.300  mm,  and  then 
to  slow  gradually  to  co/9*0.  (The  small  non-monotonic  feature  at  shallow  depths 
is  a  numerical  artifact  caused  by  applying  the  max(')  function  to  a  peak  that  is 
broad.)  The  numerically  determined  fast  speed  cq  and  the  slow  speed  co/9.0  agree 
-  quantitatively  with  analytical  results  in  Subsection  2A,  and  also  substantiate  the 
results  in  the  next-to-last  paragraph  of  Subsection  2B. 

The  last  paragraph  of  Subsection  2B  concerns  elemental  pulses  which 

are  zero  except  on  a  single  time  interval,  during  which  they  are  either  strictly  posi¬ 
tive  or  strictly  negative.  For  example,  the  Green  function  G£[z^t)  is  an  elemental 
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Fig.  6.  The  La  norm  and  the  Loo  norm  Fi{z)  of  the  Green  function  Gb,45»  (z,  0* 

The  paragraph  below  (3.13)  relates  these  norma  to  upper  bounds  for  obliquely 
incident  pnlaea. 


Time  (psec) 


Depth  (meter) 


8.  The  Lj  norm  F^^z)  and  the  Loo  norm  Fs{z)  of  the  Green  function  Gh,4s*  (2,  •) 
for  magnetic  fields  produced  by  pulses  that  axe  incident  on  water  at  an  angle 
of  45®.  The  last  paragraph  of  Section  3  relates  this  figure  to  upper  bounds 
.  fotdtH:  • 


poilse  for  aU  of  the  depths  graphed  in  Fig.  1.  The  statements  in  Subsection  2B  that 
concern  the  propagation  of  elemental  pulses  aU  rely  on  the  Mowing  approxima¬ 
tion:  If  an  elemental  pulse  /brief  (0  is  much  briefer  than  an  elemental  pulse 

and  if  /brief (t)  is  concentrated  at  the  time  to,  then  l/oda/b 

|/od^/long(t  -  s)/brief(s)l  «  ll/brief (OlM/longC*  “  io)|.  That  approximation  is  asso¬ 
ciated  with  the  concept  of  ^-sequence  functions  in  the  elementary  theory  of  Dirac 
delta  functions.  We  also  note  that  the  exponential  term  in  (2.2)  quickly  becomes 
negligible  because  it  dec^s  by  a  factor  of  132  across  each  0.300  mm  interval  of 
depth.  These  three  results— the  negligible  exponential,  the  approximation  of  the 
convolution,  and  the  observation  that  Ge  is  an  elemental  pulse— together  with 
(2.2)  and  Fig.  1  yield  the  brief-pulse  result  in  Subsection  2B,  which  has  the  term 
lfcliGE{z,t)-  The  long-pulse  result  in  that  subsection  uses  the  additional  obser¬ 
vation,  below  (3.8),  that  i(?E(2,-)lli  is  approximately  constant,  and  the  time  shift 
in  the  long-pulse  result  in  Subsection  2B  also  relies  on  Fig.  7  and  results  from  the 
previous  paragraph. 

In  a  final  matter  we  would  like  to  make  it  clear  that  we  do  not  know  what 
are  the  medical  effects  of  isolated  pulses.  We  have,  however,  developed  meth¬ 
ods  that  are  flexible  enough  that  they  may  be  useful  once  the  medical  effects  are 
known.  Although  we  originaUy  developed  the  upper-bound  method  along  the  lines 
of  peak  amplitudes  and  power  densities,  following  Ref.  5,  the  method  can  eas¬ 
ily  be  extended  to,  say,  the  time  derivative  of  the  magnetic  field  if*(z,i).  To  , 
demonstrate  the  extension,  we  computed  the  Green  function  Gn  in  ifx(0,z,t)  = 

_  |^e-“V(i  -  ^/Coo)  +  -  ^)<^h(2,3)]  /(/^ocoo)  using  the  methods  of  Refs.  13- 

17,  and  we  validated  the  computation  as  described  in  Section  4.  The  Green-function 
representation  for  Hx  yields  the  following  analog  of  (3.6). 

<  e-“iat/(-)||p  +  |/(0)|1|Gh(z,-)11p  +  I1Gh(^, •)llrll5i/(-)ll„  ^here 
_  1  and  1  <  p,  q,  i*  <  oo.  The  norms  F5(z)  =  11Gh,45«  (2,  OIU 
and  Fsiz)  =  iGH.45* (2, Oioo  for  that  inequality  are  graphed  in  Fig.  8  for  the  45°- 
incidence,  magnetic  Green  function  Gh,45®‘ 
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4.  NUMERICAL  VALIDATION  AND  LABORA¬ 
TORY  EXPERIMENTS 

We  validated  the  numerical  results  obtained  with  the  Green’s  function  approach  by 
comparing  the  computed  electric  and  magnetic  fields  in  a  Debye  half-space  (Section  2, 
^inc  “  0)  to  those  computed  with  a  finite  difference  method. Figure  9  shows 
electric  field  time-traces  computed  at  three  spatial  locations  in  the  half-space  due  to 
a  Tp  =  40-psec  square  pulse  of  initial  amplitude  1  V /m.  We  note  excellent  agreement 
to  within  a  width  of  the  line  over  an  amplitude  scale  of  10  orders  of  magnitude.  This 
indicates  a  better  than  1-paxt-per-billion  agreement.  (The  small  difference  increases  to 
two-pen-strokes’  width  at  the  shallowest  depth,  but  only  after  the  field  itself  decreases 
by  a  factor  one  thousand.)  Also,  the  speeds  of  the  first  arrival  and  of  the  peak  of  the 
response  can  be  deduced  from  this  graph.  We  see  that  the  first  arrival  bccurs  with 
speed  Co  =  c,  while  the  peak  of  the  response  arrives  with  the  speed  Ci  =  0.1116c  as 
predicted  in  Section  2.  Figure  10  shows  a  comparison  of  magnetic  fields,  computed 
with  the  two  methods  at  three  depths  in  the  half-space,  for  a  square-modulated 
sinusoidal  pulse  of  Tp  =  50-psec  duration  and  carrier  frequency  80  GHz.  Again 
sinailarly  excellent  agreement  is  noted. 

Next,  we  confirm  the  analytical  results  presented  in  Subsection  2 A  by  solving, 
with  the  Green’s  function  method,  for  the  impulse  response,  /(i  ~  ~  ^(^  “ 

of  the  Debye  half-space  described  there  for  (f>inc  =  0  (^  =  t)-  To  determine  the  two 
speeds  predicted  by  (2.1)  the  peak  of  the  impulse  response  was  obtained  from  the 
simulation  results.  The  temporal  versus  the  spatial  location  of  the  peak’s  occurrence 
is  graphed  in  Figure  7.  The  slope  of  the  graph  is  the  reciprocal  of  the  speed  of  the  peak 
of  the  response.  The  relative  unimportance  of  the  characteristic  z  =  CqI  (equivalently, 
the  wavefront  speed)  with  respect  to  the  sub-characteristic  z  =  Cyt  is  immediately 
evident.  The  value  of  the  slope  is  given  in  Figure  7  for  two  ranges  of  depth.  We  see 
that  for  depths  of  O(cot)  the  response  travels  with  the  speed  Cq,  and  then  slows  down 
in  another  O(cor)  interval.  This  additional  interval  will  be  much  smaller  or  altogether 
absent  for  band-limited  pulses.  Eventually,  the  pulse  travels  with  the  speed  cy  for 
the  remcLindev  of  the  half-space.  The  diffusive  character  of  the  main  disturbance  in 
the  half-space  is  exemplified  in  Figure  11  where  the  numerically  obtained  max{E^}  is 
compared  with  the  analytically  predicted  behavior  (1/ \/2)  as  a  function  of  depth.  The 


Electric  Field  (V/m) 


Fig.  9.  The  electric  fields,  at  three  depths  in  water,  that  result  from  a  normally 
incident  40-p8ec  square  pulse  of  1-V/m  amplitude.  At  the  three  depths 
graphed,  the  fields  axe  precisely  0  until  5.46  psec,  8.20  psec,  2ind  9.60  psec 
respectively. 
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Fig.  10.  The  magnetic  fields  that  result  from  a  normally  incident  4-cycle  1-V/m- 
amplitude  square-modulated  sinusoid  of  50-psec  duration. 
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constant  of  proportionality  for  the  \|^/z  prediction  was  fit  to  the  data  at  the  center- 
point  (4  mm)  of  the  graph.  We  note  that  2dthough  the  slow  speed  is  achieved  very 
early  on  (z  ~0.1  mm),  it  takes  some  time  for  the  pulse  to  start  diffusing  [z  ~  2.5  mm). 
This  delay  depends  on  the  frequency  content  of  the  incident  signzd,  and  will  be  shorter 
for  band-limited  pulses.  From  the  discussion  and  figures  it  is  evident  that  the  envelope 
and  duration  of  the  incident  pulse  control  the  magnitude  and  nature  of  the  response 
of  the  medium.  Any  high-frequency  carrier  component  decays  exponentially  with 

depth.  Similar  behavior  was  observed  in  numerical  simulations  with  square  pulses  of 

« 

various  durations.  All  long  duration  pulses  were  found  to  travel  unattenuated  with 
amplitude  T(u;  =  0)  X  max{E'”‘},  and  square  modulated  sinusoidal  pulses  of  various 
carrier  frequencies  and  durations  (1  to  10  cycles)  behaved  like  the  Green’s  function 
with  the  Ccirrier  component  exponenticiUy  smziU. 

We  also  checked  the  sensitivity  of  our  calculations  to  the  way  in  which  we 
modeled  the  water  data  of  Ref.  2.  Reference  25  points  out  that  many  other  data  sets 
are  available  for  water,  so  we  were  satisfied  with  any  model  that  fit  the  data  in  Ref.  2 
to  within  10%.  In  peirticulcir,  we  examined  severed  Debye-model  fits  to  the  data  in 
Ref.  2  for  frequencies  that  are  below  100  GHz  and  for  which  the  data  cire  also  smd, 
in  the  reference,  to  be  reliable.  Our  experience  in  fitting  those  data  is  that  Coo  can 
be  taken  to  be  any  number  from  1.00  through  10.0,  and  then  values  for  the  two  other 
Debye  parameters,  e,  and  r,  can  be  found  that  fit  the  data  to  within  10%.  In  that 
sense,  the  value  of  Coo  is  somewhat  arbitrary;  for  most  of  our  simulations  we  chose 
£00  =  1.00,  which  is  consistent  with  an  assumption  in  Ref.  2.  The  corresponding 
Debye  parameters  eire  £„  =  1-00  and  e,  =  80.35  and  r  =  8.13  x  10”^^s“^  in  the 
notation  of  Subsection  2A,  or,  equivalently,  a  =  9.76  x  lO^^s"^  and  b  =  1.23  x  10*'s“^ 
in  the  notation  of  (3.13).  This  fit  is  referred  to  as  Model  1  in  Fig.  12;  it  is  the  fit  that 
is  used  in  most  of  the  numerical  computations  in  this  paper.  Another  model  that  fits 
the  water  data  to  within  10%  accuracy  is  defined  by  Coo  =  5.50  and  e,  =  78.20  and 
r  =  8.1  X  10~^^s~*,  or,  alternatively,  a  =  8.98  x  10^-^s“^  and  b  =  1.23  x  10^*s“^.  This 
second  model  is  used  only  in  Fig.  12,  where  it  is  called  Model  2.  The  propagation  of  an 
incident  5-cycle  8-GHz  1-V/m-amplitude  square-modulated  sinusoid  was  computed 
for  these  two  models.  Figure  12  shows  the  resulting  electric  field  produced  at  a  9.75- 
mm  depth  for  both  models.  We  compared  the  electric  fields  at  31  other  representative 
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Electric  Field  (V/m) 


depths,  from  0  cm  through  4  cm,  and  observed  similar  agreement,  commensurate  with 
two  different  10%-accuracy  fits  to  the  water  data.  A  similar  degree  of  agreement  was 
also  seen  for  the  two  models  in  computations  for  the  following  incident  fields:  (1) 
A  1-cycle  8-GHz  squzire-modulated  sinusoid;  (2)  A  1-cycle  10-GHz  squ2Lre-modulated 
sinusoid;  (3)  A  50-psec  square  pulse;  and  (4)  A  100-psec  square  pulse.  These  results 
are  numerical  evidence  that  the  computed  internal  fields  cire  stable  as  the  water  model 
is  perturbed  by  about  10%. 

Now  we  wish  to  provide  an  explanation,  based  on  the  understanding  devel¬ 
oped  herein,  of  the  observations  in  Ref.  3  whereas  no  significant  differences  in  SAR 
(Specific  Absorption  Rate)  distribution  between  pulsed  and  CW  (Continuous  Wave) 
exposures  were  measured  for  a  MEM  (Muscle  Equivalent  Material)  exhibiting  two  re- 
Izocation  times.  The  dielectric  model  was  composed  of  two  Debye  mechanisms,  which 
fit  experimentally  determined  permittivity  and  conductivity  data  for  MEM  at  2.07, 
2.8,  5.6,  and  9.3  GHz.  The  two  relaxation  times  were  n  =  6.63  psec,  and  T2  =  83.7 
psec.  The  material  was  illuminated  with  a  train  of  squaxe-modulated  pulses  of  var¬ 
ious  durations  and  repetition  rates.  We  will  explain  the  observations  in  Ref.  3  for 
the  smallest  pulse  repetition  rate  used  (200  pulses  per  second)  for  which  the  pulse 
duration  was  0.5  /xsec.  All  other  results  in  Ref.  3  with  different  pulse  settings  are 
similsurly  explained.  For  this  incident  signal  the  author  of  Ref.  3  used  a  carrier  of 
5,6  GHz.  Thus,  the  pulse  duration  was  Tp  =  5  x  10”^  sec,  the  quiescent  interval 
between  pulses  was  Tq  =  5  X  10“^  sec,  and  each  pulse  contained  2800  cycles  of  car¬ 
rier.  Further,  it  happened  that  Tp  >>  max{r}  =  r2,  thus  the  medium  would  not 
respond  in  a  dispersive  manner,  rather  it  exhibited  an  effective  relative  permittivity 
of  e]  +  e]  --  600  =  42.4  +  15.8  —  4.3  =  53.9  to  the  envelope  of  the  pulses.  With  these 
parameter  values  in  mind  we  expect  the  medium  to  sense  a  CW  signal  of  carrier  5.6 
GHz,  even  in  the  pulsed  case.  Since  the  MEM  has  low  heat  diffusivity,  i.e.,  it  takes 
about  40  sec  for  a  temperature  change  of  0.04  to  occur,  its  temperature  will  change 
immeasurably  in  the  time  Tq  between  pulses.  Consequently,  the  temperature  will  re¬ 
main  constant  until  the  next  pulse  arrives  again  to  be  seen  by  the  medium  as  part 
of  a  CW  exposure,  hence  the  observed  indistinguishability  of  the  CW  SAR  versus 
the  pulsed  SAR.  As  the  carrier  frequency  is  increased  the  CW  vs.  pulsed  exposure 
SARs  should  start  disagreeing  at  smaller  depths.  This  is  related  to  the  frequency- 
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dependent  skin  depth  for  the  carrier  component  which  is  9.7  mm  at  5.6  GHz.  AU 
the  measurements  of  SAR  were  obtained  at  depths  weU  within  the  skin  depth.  The 
effects  of  the  pulsing  should  be  observable  at  depth  greater  than  CoTi  ~  1.2  cm  since 
then  the  carrier  component  will  have  decayed  sufficiently  (it  is  exponentially  decaying 
as  in  the  CW  case)  so  the  remaining  field  wUl  be  due  to  the  square  envelope  and  will 

behave  diffusively. 

The  results  presented  in  our  paper  also  help  in  accurately  predetermining  sam- 
pie  thickness  to  be  used  in  single  and  total  transmission  Time  Domain  Spectroscopy 
(TDS)  studies  such  as  those  presented  in  Ref.  26.  In  the  single  transmission  approach 
one  studies  the  first  arrival  through  a  long  sample  so  that  the  highest  frequency 
components  will  have  decayed  sufficiently  in  order  not  to  mask  the  lower  frequency 
components  which,  as  we  showed  in  Section  2,  are  significantly  slower  (for  water 
Cl  ~  Co/9).  The  first  suggestion  arising  from  the  analysis  of  Section  2  (verified  by 
the  numerical  simulations)  is  that  the  sample  can  be  as  short  as  2c<,r  when  one  is 
interested  in  measuring  the  static  permittivity  of  a  Debye-like  material  with  single 
transmission  TDS.  On  the  other  hand,  in  the  total  transmission  approach  the  first 
arrival  time  through  a  short  sample  is  used  to  best  measure  the  infinite-frequency 
permittivity  of  the  material  under  test.  Therefore,  a  sample  length  shorter  than  c^t 
should  be  appropriate  in  order  to  capture  the  time  of  arrival  of  the  highest  frequencies 
and  thus  determine  (The  next-to-last  paragraph  of  Section  2  has  related  com¬ 
ments.)  For  the  test  case  presented  in  the  results  section  of  Ref.  26  [Eqn.  (8)  there] 
it  happens  that  e,  =  17.3,  £„  =  3.3,  r  =  460  psec.  Thus,  the  highest  frequency 
components  travel  with  a  speed  Co  =  0.5505c,  the  lowest  frequency  components  travel 
with  speed  Ci  =  0.2673c,  and  c^r  =  7.6  cm.  The  experimenters  could  have  used  a 
slab  of  material  of  at  least  15  cm  to  reliably  measure  e.,  and  a  slab  thickness  of  at 
most  7  cm  to  measure  Coo,  instead  of  using  50  cm  and  2.7  cm  respectively  to  do  the 
measurements.  In  Figure  13  we  show  the  x-t  location  of  the  peak  of  the  impulse 
response  for  this  medium.  (The  two  data  points  in  one  corner  of  the  graph  are  minor 
numerical  cirtifacts  that  are  related  to  the  similar  numerical  artifacts  in  Fig.  7.)  The 
graph  as  a  whole  confirms  our  predictions  of  the  sample  lengths  with  which  TDS  will 
be  most  successful  in  measuring  the  permittivities. 
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Fig.  13.  The  time  of  arrival  of  the  peak  of  the  Green  function  Gsiz.t)  at  various 
depths  in  the.  medium  of  Ref.  26. 
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5.  CONCLUSION 

There  are  many  generalized  methods  for  computing  the  response  of  any  single  dis¬ 
persive  medium  to  any  single  incident  pulse.  This  paper  has  contributed  nothing 
along  those  lines;  instead,  we  have  developed  several  related  results  that  are  general 
in  a  different  way.  We  have  found,  for  instance,  that  the  power  density  and  the 
peak  amplitude  of  an  incident  pulse  place  upper  bounds  on  the  peak  amplitude  suf¬ 
fered  inside  a  dispersive  medium — independent  of  the  other  details  of  the  incident 
pulse.  Our  Green-function  results  reduce  the  computation  of  such  upper  bounds 
to  little  more  than  a  calciilator  exercise,  and  we  have  given  numerical  examples  in 
which  these  sharp  upper  bounds  are  almost  attained.  We  have  reported  similar 
upper-bound  estimates  for  a  quantity  related  to  power  density,  and  for  the  time 
derivative  of  the  magnetic  field.  That  time  derivative,  dtH,  is  largely  responsible 
for  electromotive-force  currents  in  circuit-like  structures;  it  is  especially  large  for 
short-risetime  pulses.  Such  upper  bounds  could,  potentially,  help  in  the  regulation 
of  electromagnetic  interference  or  damage  produced  in  dispersive  media. 

Although  our  methods  apply  to  dispersive  media  generally,  we  have  used 
a  Debye  model  for  microwave-pulse  propagation  in  water  as  a  specific  numerical 
example.  For  that  water  example,  we  reported  a  factor-of-nine  effect  in  the  wave 
speeds  that  seems  to  have  been  unnoticed  until  now,  and  we  explained  this  large 
effect  analytically  using  PDEs.  The  PDE  analysis  also  yielded  simple  short-pulse 
and  long-pulse  approximations,  as  did  an  analysis  of  the  numerical  Green  functions 
involved  in  the  upper-bound  concepts  described  earlier.  We  studied  rates  of  decay 
in  Subsection  2A,  and  approximations  for  detailed  time  traces  are  given  in  Subsec¬ 
tion  2B  and  Eqns.  (3.2)  and  (3.4).  Section  4  uses  that  work  to  explain  the  results 
of  some  laboratory  experiments,  and  to  offer  suggestions  for  future  experiments. 
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APPENDIX  A:  ERROR  ESTIMATES 

This  appendix  develops  an  easily  used  method  for  estimating  the  percentage  error 
in  grid-dependent  computations.  The  estimates  are  validated  numerically  in  some 
cases  for  which  exact  solutions  are  known,  and  the  estimates  are  also  used  in  cases 
for  which  exact  solutions  are  not  known.  These  estimates  use  only  one  of  the  many 
definitions  of  relative  error  for  functions  of  two  variables,  but  the  estimates  can  easily 
bei  adapted  to  other  measures  of  relative  error  and  to  functions  of  more  than  two 
variables.  We  were  motivated  to  estimate  percentage  errors  because  quantitative 
estimates  of  uncertainty  could  help  in  setting  safety  standards  or  in  the  use  of  other 
computations  for  which  measures  of  numerical  uncertainty  are  vital. 

The  error  estimates  in  this  paper  are  empirical.  The  estimates  are  obtained 
by  computing  with  severzd  different  grid  sizes,  noticing  a  pattern  in  the  relative 
errors  of  the  different  computer  runs,  and  using  that  inferred  pattern  to  sum  the 
series  on  the  right  side  of  the  inequality  ||/c  —  /c  i  oo  ^ 

/a loo  +  The  concepts  of  rate  and  order  of  convergence  are  not  used  in  these 
empirical  estimates.  The  logic  behind  the  estimates  is  slightly  intricate,  but  the 
resulting  method  uses  only  least  upper  bounds  and  geometric  means  and  geometric 
series,  which  are  extraordinarily  simple  concepts. 

We  will  estimate  relative  errors  for  functions  of  two  variables.  The  error  in 
a  computed  solution  /c(-2)t)  relative  to  the  exact  solution  /e(z,t)  can  be  defined  as 


i/e  (^,-)  loo 


{Al) 


The  relative  error  (Al)  is  a  function  of  the  depth  z.  For  instance,  if  the  relative 
error  is  known  to  be  less  than  1%  over  a  range  of  depths  z,  then,  for  those  depths, 
the  exact  solution  /e(z,t)  can  differ  from  the  computed  solution  fc{z,t)  by  no  more 
than  1%  of  the  peak  value  (with  respect  to  t)  of  |/e(z,t)|.  The  just-mentioned  peak 
value  of  |/e(zj")|  is  unknown  in  zdl  practical  cases  in  which  the  exact  solution  /e  is 
unknown;  however,  in  the  present  hypothetical  case  of  1%  relative  error,  the  peak 
value  of  the  unknown  quantity  |/e(^j  *)!  differ  by  no  more  than  1%  from  the  peak 
value  of  the  known  quantity  |/c{^}  *)1'  ^  brief  derivation  involving  a  geometric  series 
then  shows  that,  in  the  case  of  1%  relative  error,  the  graph  of  the  exact  solution 
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is  guaranteed  to  fall  between  the  graph  of  fc{z,t)  +  l-01%J/c('Z, •)||oo  and 
the  graph  of /c(i,t)-1.0T%||/c(z,-)Ioo.  The  quantities  ±1.01%|/c(2,  OIU  are  caUed 
error  bars. 

We  will  now  go  step  by  step  through  the  error  estimates  illustrated  in  Fig.  14. 
Those  estimates  are  for  Green-function  computations  of  waves  normally  incident  on 
water.  The  computations  were  done  for  five  different  grids,  and  the  results  are 
called  in  order  of  the  increasing  fineness  of  the  grids. 

Because  of  computer  limitations,  the  finer-grid  computations  were  done  for  shal¬ 
lower  depths  than  were  the  coarser-grid  computations.  For  zero  depth,  the  light 
circles  on  the  z  =  0  axis  of  Fig.  14  represent  l/i(0,-)  -  /zCO, •)loo/|/5(0, Olloo, 
1/2(0,-)  -  /3(0,-)ico/i/5(0,-)ioo,  1/3(0,-)  -  /4(0,-)lco/l/5(0,-)ic^,  and  1/4(0,-)  - 
/5(0,-)loo/l/5(0,-)ioo,  running  from  top  to  bottom  along  the  vertical  a^s.  For  in¬ 
stance,  the  error  in  /4(0,t)  relative  to  the  most  finely  computed  result  /5(0,t)  is 
0.157%.  The  even  spacing  of  the  four  light  circles  along  the  logarithmic,  vertical 
axis  of  Fig.  14  suggests  that  the  four  relative  errors  are  in  geometric  progression; 
in  fact,  the  ratio  of  the  second-largest  relative  error  to  the  largest  relative  error 
is  0.368,  the  ratio  of  the  next  two  smaiUer  relative  errors  is  0.284,  the  ratio  of  the 
two  smallest  relative  errors  is  0.268,  and  the  geometric  mean  of  the  three  previous 
numbers  is  0.304.  For  those  reasons  we  assume 

l/n-n(0,-)-/n+2(0,-)i~  _  (0.304)^^  — 

l/s(0,-)loo  1/5(0, -)ico 

even  for  hypothetical  results,  such  as  /loo,  which  would  come  from  computations 
involving  a  mucli-finer  grid  than  was  used  for  the  actual  computation  of  /s.  We 
also  assume  that  the  exact  result  is  the  n—>cx)  limit  of The  triangle 

inequality  for  norms  then  implies  ||/i  —  /« loo  ^  l/i  /2  llcxa  +  II/2  ^  /a  Icxs  +  i/s  “  A  l|oo  + 
•  •  •.  We  use  (A2)  in  the  right  side  of  the  previous  inequality  to  get  a  geometric  series, 
which  sums  to 


1/1  (Q,-)  ~  /e(0r)loo  <  ^  _ /2(0,  Oioo  _  g  (^3) 

1/5(0, -)lco  1  -  0.304  1/5(0, -)lco 
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Relative  Error  (1  =  100%) 


I 


Fig.  14.  Relative  errors  for  the  computation  of  Ge  for  water.  The  bold  marks  rep¬ 
resent  the  final  error  estimates  and  their  numerical  validations.  The  lighter 
marks  represent  intermediate  steps  in  the  computation  of  error  estimates. 
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at  z  =  0.  That  8.08%  estimated  relative  error  is  plotted  as  the  top-most  heavy  circle 
on  the  vertical  axis  of  Fig.  14.  The  analysis  was  repeated  for  the  relative  errors  of 
other  computations  yielding 


n/i(or)-/«(o.-)ic 
1/5(0,  Olco 


1/5(0,  Oico 


foj  t  =  1,2,3,4.  Those  estimated  relative  errors  are  also  plotted  as  heavy  circles  on 
the  vertical  axis  of  Fig.  14. 

Reference  9  has  an  exact,  closed  form,  modified-Bessel-function  expression 
of  the  zero-depth  Green  function  for  a  half-space  of  Debye  medium.  That  exact 
solution  /e  was  used  to  compute  the  true  relative  errors  E,ei[/i,/e]  at  z  =  0.  Those 
true  relative  errors  are  plotted  as  the  bold  Xs  on  Fig.  14.  The  true  relative  errors 
(bold  Xs)  validate  the  error  estimates  (bold  circles)  because  the  two  sets  of  errors 
match  closely;  for  instance,  the  error  estimate  (A4)  for  /3(0,-)  was  that  there  would 
be  0.844%  error  relative  to  /«,  and  the  actual  error  in  /3(0,-)  relative  to  /*  was 
0.797%.  Our  percentage-error  method  involving  geometric  means  and  geometric 
series  is  thereby  validated. 

The  relative  errors  l/i(^,-)  —  fi+i{^i')loo/lh{^i‘)loo  were  then  computed 
at  the  depths  z  =  0.480  mm,  1.04  mm,  1.52  mm,  and  2.00  mm,  for  i  =  1,2, 3,4. 
Those  results  are  plotted  as  some  of  the  light  circles  to  the  right  of  the  vertical 
axis  in  Fig.  14.  Error  estimates  for  fz{zr)  were  computed  using  the  geometric- 
series  technique  in  (A2)-(A4),  with  the  results  plotted  as  several  bold  asterisks  on 
Fig.  14.  The  ratios  of  relative  errors  and  the  geometric  means  of  the  ratios  were 
computed  separately  for  each  depth.  In  particular,  the  geometric  means  for  the  four 
depths  mentioned  in  this  paragraph  are  0.302,  0.328,  0.315,  and  0.333  for  0.480  mm 
through  2.00  mm,  consecutively.  We  will  now  explain  how  the  relative  errors  for 

depths  beyond  2.00  mm  were  computed. 

The  finest-discretization  run  was  not  computed  beyond  2.00  mm  for  reasons 
related  to  computer  resources.  Thus,  only  /i,  /z,  /s,  a-^d  fi  were  analyzed  for  the 
depths  from  2.48  mm  through  4.00  mm.  These  analyses  were  done  as  described  in 
the  previous  paragraph,  and  the  geometric  means  were  also  computed  separately 
for  each  depth.  Similarly,  the  runs  /i,  /z,  and  /z  were  analyzed  for  depths  from 
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5.04  TTim  through  7.04  mm.  In  these  computations,  and  in  all  other  computations 
described  in  this  appendix,  the  Loo  norms  were  computed  using  90  evenly  spaced 
time  points. 

The  final  result  in  Fig.  14  is  symbolized  by  the  bold  asterisks  there.  The 
result  is  that  has  errors  relative  to  /.  that  are  expected  to  fall  monotonically 

from  about  1.70%  at  z  =0.480  mm  to  about  0.799%  at  z  =4.00  mm.  The  error  in 
/3(0,i)  relative  to  /e(0,t)  is  actually  0.797%.  Figure  14  shows  that  the  relative 
errors  are  non-monotonic  from  0.00  mm  through  0.500  mm;  therefore,  we  make  no 
further  inferences  about  the  relative  errors  in  that  first  half- millimeter. 
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